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SUMMARY 


A  mathematical  model  to  predict  the  dynamics  of  a 
flexible  orbiting  platform  is  developed.  The  platform  is 
idealized  as  a  large  thin  homogeneous  square  plate,  made  up 
of  a  continuous  distribution  of  mass  points  in  a  plane.  By 
considering  the  internal  and  external  forces  acting  on  each 
generic  mass  point,  the  equations  for  the  rigid  body 
motions,  as  well  as  the  elastic  degrees  of  freedom  are 
developed.  It  assumed  that  the  elastic  motion  is  limited  to 
small  amplitudes  and  that  the  center  of  mass  follows  a 
circular  orbit.  For  small  amplitude  flexural  motion,  the 
rigid  body  and  elastic  modes  are  modeled  to  the  first  order, 
thus  linearizing  the  equations  of  motion  for  control  law 
synthesis. 

Under  the  assumption  that  the  linear  system  is 
completely  observable,  the  optimal  control  laws  are 
developed  for  the  case  where  the  observational  data  is 
collected  on  a  sampled  basis,  i.e.,  a  discrete  time  data 
system. 

The  attitude  and  shape  control  can  be  achieved  by 
placing  point  thrust  actuators  perpendicular  to  the  main 
surface  and  the  edge  of  the  plate.  Their  effects  on  the 
motion  are  modeled  to  the  first  order.  Controllability  for 
the  system  is  verified  for  two  sets  of  actuator  locations. 

An  application  of  the  linear  quadratic  regulator  technique 
in  a  discrete-time  domain  yields  the  optimum  control  law 
feedback  gains. 

A  comparison  of  the  performance  of  the  different  sets 
of  actuator  locations  results  in  the  best  choice  of  actuator 
positioning.  Parametric  studies  are  conducted  to  show  the 
effect  of  varying  the  state  penalty  matrix,  control  penalty 
matrix,  and  the  sampling  period  on  the  transient  performance 
of  the  system. 
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CHAPTER  1  -  INTRODUCTION 


The  Solar  System  is  our  extended  home.  Through  the  U5C 
of  space  technology  mankind's  expansion  outward  from  Earth 
to  other  worlds  becomes  feasible.  The  opportunity  to 
stimulate  individual  initiative  and  free  enterprise  in  space- 
will  surface  through  the  development  of  new  lands. 
Historically,  when  the  power  of  the  human  intellect  combines 
with  abundant  energy  and  rich  material  resources  wealth  is 
created.  On  the  space  frontier,  new  wealth  can  be  created  to 
benefit  the  entire  human  community  by  combining  the  energy 
of  the  Sun  with  materials  left  in  space  during  the 
formulation  of  the  Solar  System.  Mankind's  reach  will  extend 
in  science,  industry,  and  the  settlement  of  space  with  the 
correct  combinations  of:  vigor  and  continuity,  the  elements 
of  scientific  research,  technological  advances,  the 
discovery  and  development  of  new  resources  in  space  and  the 
provisions  of  essential  industries  and  systems  .  Government 
investments  will  generate,  in  value,  financial  returns  many 
times  its  initial  cost  to  the  benefit  of  all. 

To  meet  the  challenge  of  the  space  frontier,  the 
National  Commission  on  Space,  has  proposed  a  step  by  step 
program  to  open  the  inner  Solar  System  for:  exploration, 
basic  and  applied  research,  resource  development,  and  human 
operations . With  the  advent  of  the  Space  Shuttle  as  a 
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reliable,  affordable  transportation  system,  the  preliminary 
steps  of  acquiring  a  network  of  outposts  in  space  can  be 
undertaken.  By  following  a  systematic  program  -with  minimum 
risk  and  funding-  a  progressive  path  for  future  space 
activities  can  occur.  This  program's  structure  will  be  in 
accordance  with  the  inner  Solar  System's  natural 
characteristics:  energy,  distance,  signal  delay  time,  and 
availability  of  resources. 

An  outpost  actually  consists  of  one  or  several  large 
platforms  with  connecting  appendages,  i.e.,  solar  panels, 
radar  dishes,  habitation  modules,  skylabs,  etc..  There  are 
major  differences  between  satellites  and  large  platforms; 
size  and  capabilities  are  two  such  differences.  Large 
platforms  are  well  defined  by  Cuneo  and  Williams,  '  as  a 
system  which  provides  basic  services  to  a  changing  set  of 
activities.  '  The  capability  of  service  -either  updating 
payloads,  or  performing  repairs,  or  replacing  degraded 
modules,  and/or  replacing  consumables-  is  probably  the  most 
common  aspect  attributed  to  platforms.  The  primary  purposes 
of  a  platform  is  to  provide  shared  support  for  multiple 
payloads  and  to  provide  connectivity.  The  reasons  for 
platforms,  are  to  obtain:  (1)  the  economies  of  scale  which 
come  from  shared  support,  and  (2)  the  new  and  improved 
services  which  come  from  connect ivity . 


The  forementioned  network  of  outposts  would  have 
several  different  locations:  low  earth  orbit  (LEO) ; 
geostationary  orbit  (GEO);  lunar  (surface  and/or  orbit); 

Mars  and  its  asteroids.  The  outposts  of  interest  for  this 
thesis  will  be  located  in  LEO,  approximately  250  nautical 
miles  from  the  Earth's  surface.  Low  earth  orbits  are  those 
just  beyond  the  Earth's  atmosphere  and  are  the  easiest  to 
reach  from  Earth.  This  orbit  provides  both  a  close  proximity 
orbital  view  of  Earth  and  a  window  for  observation  of  the 
Universe.  Freedom  from  strong  gravitational  effects  allows 
experiments  that  would  be  impossible  to  conduct  on  Earth's 
surface  and  facilitates  the  construction  of  large  structures 
of  low  mass.  Earth  provides  a  sheltering  skirt  of  magnetic 
field  that  protects  us  from  solar  flare  radiation. 

Planetary  landings  are  costly  in  terms  of  propellant 
requirements  for  the  descent,  but  the  access  of  surface 
materials  becomes  an  invaluable  resource.  When  lifting 
payloads  into  orbit,  away  from  Earth's  gravitational  field, 
we  expend  energy;  to  overcome  the  Earth's  gravitational 
clutch,  the  rockets  must  attain  speeds  to  lift  a  payload 
free  of  Earth's  pull.  We  must  expend  the  same  amount  of 
energy  necessary  to  haul  that  same  payload  influenced  by  the 
full  force  of  gravity  to  a  height  of  4,000  miles.  To  reach  a 
nearer  goal  of  low  Earth  orbit  -where  rockets  and  their 
payloads  achieve  a  balancing  act,  while  skimming  above 
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Earth’s  atmosphere-  we  must  spend  about  half  as  much  energy, 
equivalent  to  climbing  a  mountain  2,000  miles  high.*'^  Once 
in  'free  space',  away  from  planets  and  moons,  large 
distances  can  be  traveled  with  modest  expenditures  of 
energy.  Other  gains  from  working  in  an  orbiting  space 
environment  include:  full-time  solar  energy,  which  is 
valuable  for  industrial  processing,  and  microgravity,  which 
is  advantageous  for  building  large  space  structures,  i.e. 
Space  Station,  Variable-G  Research  Facility,  and  a  Mars 
Interplanetary  Vehicle.  Once  these  facilities  are  built, 
research  data  can  be  accumulated  to  understand  how  the 
absence  of  gravity  affects  the  fabrication  of  faultless 
materials,  physical  conditions,  and  motor  skills  of 
humans. 

Early  industrial  production  in  space  may  be  best 
achieved  by  transporting  raw  materials  from  the  Moon  or 
Earth  to  orbiting  platforms  processing  and  fabricating 
finished  products  via  robotic  factories  powered  by 
continuous  solar  energy. 

The  first  space  enterprise  to  reach  economic  viability 
was  satellite  communications.  Once  in  orbit,  communication 
satellites  lock  into  geostationary  position  and  relay 
electronic  messages,  telephone  calls,  electronic  mail,  and 
television  broadcasts.  Future  developments  in  space-based 


communications  and  information  systems  will  revolut ioni Ee 
our  daily  lives.'*’  The  proposed  deployment  of  large  plate  or 
dish  shape  orbiting  structures  will  make  it  possible  to 
equip  a  car,  boat,  or  airplane  with  a  receiver  and  a  display 
to  pinpoint  its  exact  location  by  satellite,  allowing  the 
provision  of  navigation,  collision  warning,  fleet  dispatch, 
emergency  location,  and  two-way  communication  via  satellite. 
These  seirvices  can  even  be  provided  to  small  hand-held 
terminals  powered  by  penlight  batteries. 

Currently  in  early  stages  of  development  is  the  remote 
sensing  from  low  Earth  orbiting  satellites  and/or 
structures.  From  the  vantage  point  of  space  they  enhance  the 
ability  to  observe  and  produce  specialized  maps,  that  help 
facilitate  the  management  of  crops  and  mineral  resources  as 
well  as  forecast  potentially  destructive  phenomena  to 
forests,  fisheries,  pollution,  and  water  resources. 

One  highly  prospective  space  enterprise  would  -if 
technically  and  economically  feasible-  satisfy  all  of  the 
ideal  conditions.  This  enterprise  would  provide  energy  for 
Earth  from  orbiting  structures  that  are  intercepting  solar 
energy.'^’  The  Soviet  Union  has  already  announced  the  goal  of 
building  the  first  solar  powered  satellite  to  supply  energy 
to  the  Earth  in  the  1990 's.  To  capture  such  a  market  would 
have  substantial  impact  on  the  world's  energy  problem."’ 


Using  the  capabilities  of  the  Space  Shuttle  many  future 
missions  have  been  proposed  based  on  the  deployment  of  large 
space  structures.  Most  of  the  satellites  that  have  been 
launched  so  far  consist  of  a  massive  central  rigid  body, 
with  characteristic  dimensions  of  the  order  of  a  few  meters, 
attached  to  light  rigid  or  flexible  appendages  usually 
characterized  by  dimensions  of  not  more  than  tens  of  meters. 
The  natural  frequencies  associated  with  such  flexible  parts 
are  normally  several  orders  of  magnitude  greater  than  the 
orbital  frequency  and  the  frequencies  of  the  rigid  body 
rotational  motions.  In  contrast  to  the  existing  satellite 
systems,  the  proposed  large  space  systems'  entire  structure 
will  be  considered  flexible.  With  the  inherent  size  and 
necessary  low  weight  to  area  ratio,  the  structural 
frequencies  in  the  range  of  1/100  Hertz  or  less  may  be 
considered  in  the  study  of  the  dynamics  and  control  of  these 
systems.  For  these  proposed  missions  the  operational 
considerations  define  stringent  accuracy  limits  (possibly  of 
the  order  of  millimeters,  for  a  typical  structure  of  lOO  m) 
on  the  shape  control  of  these  structures.^®’  To  satisfy  these 
requirements  and  others,  both  shape  and  orientation  of  the 
orbiting  system  should  be  controllable. 

Often  the  optimal  control  laws  for  these  future  systems 
are  developed  under  the  assumption  that  the  state  vector  is 
observed  directly  or  the  state  information  can  be  estimated 


on  a  continuous  basis.  However,  for  future  applications  the 
observational  data  will  often  be  collected  on  a  sampled 
basis,  creating  a  discrete  tine  data  system.  The  amount  of 
information  collected  may  be  reduced  and  the  format  of  data 
input  may  be  acquired  more  conveniently.  The  case  presented 
here  will  have  the  characteristics  of  being  completely 
observable,  with  an  addressed  deterministic  system  (i.e.,  no 
random  noise  nor  sensor  system  dynamics  will  be  considered). 

It  will  be  useful,  and  timely,  to  study  the  control 
problem  of  large  flexible  orbiting  space  structural  systems 
with  discrete-time  observational  data.  The  development  of 
modem  control  theory  and  technology  provides  a  strong  tool 
for  solving  this  kind  of  engineering  problem.  The  LQR 
regulator  technique  is  that  strong  tool  for  synthesizing 
linear  system  control  laws.^’^  The  LQR  strategy  can  provide 
acceptable  control  performance  once  the  state  and  control 
penalty  matrices  are  properly  selected.  It  does  not  restrict 
the  number  of  actuators  to  be  equal  to  the  number  of  degrees 
of  freedom  in  the  system.  Although,  the  LQR  method  has  been 
developed  and  widely  applied,  it  is  still  not  an  easy  task 
to  apply  it  to  the  engineering  design  for  the  control  of 
large  space  structural  systems,  especially  for  systems  with 
sampled  data  input. There  are  still  many  specific 
problems  to  be  investigated.  This  present  work  represents  an 
initial  effort  toward  understanding  the  dynamics  and  control 


of  large  flexible  structures  in  a  discrete  time  domain. 

By  using  a  continuum  approach,  Santini,  developed  a 
mathematical  formulation  for  predicting  the  motion  of  a 
general  orbiting  flexible  body.''”’  The  formulation  is  based 
on  the  following  assumptions  that;  (1)  the  elastic 
deformations  are  in  the  linear  ranges  so  that  the 
displacements  can  be  expressed  as  superpositions  of  the 
various  modes;  and  (2)  the  linear  characteristic  dimensions 
of  the  structure  are  assumed  much  smal  er  than  the  orbital 
radius.  The  effects  of  higher  harmonics  in  the  Earth's 
gravitational  potential  are  included.  However,  the 
formulation  has  a  slight  drawback  in  that  the  elastic  modal 
shapes  are  expressed  only  in  terms  of  Cartesian  components. 
This  can  be  remedied  by  redeveloping  the  formulation  using 
vector  calculus,  which  was  done  by  Kumar.'®’ 

For  simple  structures,  such  as  beams,  plates,  and 
shallow  spheres  moving  in  orbit,  modified  versions  of 
Santini ' s  formulation  allow  the  development  of  the 
translational,  rotational,  and  elastic  equations  of  motion. 
In  order  to  gain  insight  into  the  dynamics  and  control  of 
the  proposed  large  flexible  platform  system,  the  formulation 
and  manipulation  of  the  equations  of  motion  for  a  free-free 
beam  were  studied.””^'”'  Cf  particular  interest  was  a  beam 
which,  in  equilibrium,  has  its  longitudinal  axis  aligned 


along  the  vertical. 


Assuming  the  center  of  mass  follows  a  circular  orbit 
and  the  pitch  and  the  flexural  deformations  occur  only 
within  the  orbital  plane,  it  is  seen  that  the  pitch  motion 
does  not  influence  the  elastic  motion.  The  pitch  and  the 
elastic  modes  are  decoupled  for  large  values  of  the  square 
of  the  ratio  of  the  structural  modal  frequency  to  the 
orbital  angular  rate.  For  small  values  of  this  ratio  the 
elastic  motion  is  governed  by  Hill's  3-term  equation  w'hich 
can  be  approximated  by  a  Mathieu  equation.  Using  a  Mathieu 
stability  chart,  the  resulting  stability  was  considered . 

For  small  amplitude  flexural  motion,  the  rigid  body  and 
elastic  modes  were  modeled  to  the  first  order,  thus 
linearizing  the  equations  of  motion.”^’  A  parametric 
analysis  of  the  controllability  of  the  motion  of  the  beam 
about  a  nominal  orientation  for  a  discretized  system  was 
reviewed.  Extensions  were  made  to  the  LQR  formulation  and 
were  applied  to  a  thin  flexible  orbiting  plate  in  a 
continuous  time  domain. 

Many  large  space  structures  proposed  for  future  space 
applications  can  be  approximated  by  the  basic  structural 
forms  of  a  thin  plate  or  a  shallow  spherical  shell.  The 
ability  to  accurately  determine  the  frequencies  and  mode 
shapes  is  essential  for  the  analysis  and  control  of  large 


orbiting  structures.  Methods  of  describing  free  and  forced 
vibrations  of  plates  and  shallow  shells  have  been  formulated 
by  many  investigators,  A  comparative  study  of  several 
different  methods  was  reviewed  for  a  free-free  aluminum 
square  plate. 

An  extension  to  the  comparative  study  of  reference  19 
has  been  accomplished  in  this  thesis  to  include  an 
additional  technique,  GTSTRUDL.  GTSTRUDL  is  a  structural 
design  language  that  incorporates  the  finite  element  method. 
With  this  additional  technique  the  final  product  is  now  a 
comparison  of  four  different  frequency  and  mode  shape 
approximation  methods. The  methods  compared  were; 

1)  the  approximate  frequencies  and  mode  shapes  of  a 
rectangular  plate  derived  from  the  formulation  by 
Warburton  ^^2,23] . 

2)  the  analytical  results  for  a  square  plate  were 
calculated  from  a  method  by  Lemke'^^^; 

3)  the  frequencies  and  mode  shapes  were  computed 
using  a  finite  element  program,  STRUDL,  written  at 
M.I.T.  ;  and 

4}  the  frequencies  and  mode  shapes  were  computed 
using  GTSTRUDL,  an  update  version  of  STRUDL 
written  at  Georgia  Tech  . 

It  was  found  that  GTSTRUDL  obtained  better  results  than 
STRUDL  and  produced  accurate  results  for  specific  finite 


element  (see  Table  1)  input  grid  point  (node)  locations, 
whereas,  the  Warburton  and  Lemke  methods  could  only  afford 
approximate  answers. 

Attitude  and  shape  control  are  assumed  to  be  realized 
in  this  investigation  by  placing  point  thrust  actuators  in  a 
symmetrical  pattern  perpendicular  to  the  main  surface  and 
the  edge  of  the  plate.  The  placement  of  the  actuators  on 
the  main  surface  help  control  the  shape  deformation  and  the 
torgue  about  two  of  the  principal  axes.  The  placement  of  the 
actuators  along  the  edge  of  the  plate  help  control  the 
torque  about  the  third  axis.  Care  is  taken  not  to  place  the 
actuators  on  the  nodal  lines  associated  with  the  first  and 
second  transverse  vibrational  modes  or  on  the  nodal  circle 
associated  with  the  third  mode.  The  actuator's  effects  on 
the  rigid  body  and  elastic  modes  are  modeled  to  the  first 
order.  In  this  thesis  it  will  be  assumed  that  the  system  is 
completely  observable. 

The  control  laws  for  the  system  will  be  applied  to 
obtain  the  optimal  control  feedback  gains  based  on  an 
application  of  the  linear  regulator  problem  for  a  discrete 
time-data  system.  The  implementation  of  the  LQR 

procedure  will  be  accomplished  by  using  the  ORACLS  software 
routines . 


A  system  may  be  completely  controllable  in  a 
continuous-time  domain,  but  once  it  is  discretized, 
controllability  is  not  guaranteed.  To  insure  the 
controllability  of  a  discretized  data  system,  the  sampling 
period  theorem  will  be  applied  for  proper  selection  of  the 
sampling  period.  The  theorem  mandates  the  eigenvalues  of  the 
closed  loop  system  must  satisfy  certain  conditions . For 
signal  reconstruction,  the  sampling  period,  AT,  snould  be  as 
small  as  possible,  but  if  the  sampling  time  is  too  small, 
the  computational  requirements  may  exceed  the  computer 
speed . 

Under  noirmai  operation,  the  onboard  computer  estimation 
and  control  must  finish  processing  all  the  input  data  during 
one  sampling  period,  AT,  i.e.,  the  prediction  of  the  state 
variables  which  will  be  used  for  the  controller  must  be 
available  before  the  beginning  of  the  next  sampling 
sequence.  Thus  the  sampling  period  should  be  more  than  the 
minimum  computational  time  required  by  the  onboard 
microcomputer  for  the  simulation  of  each  step  in  the 
estimation  and  control  process.  The  choice  of  sampling  time 
is  also  constrained  by  the  performance  of  the  transient 
response,  i.e.,  oversh.jOt  characteristics,  settling  time, 
steady  state  RMS  errors. 


CHAPTER  2  -  DEVELOPMENT  OF  EQUATIONS  OF  MOTION  FOR  A 
FLEXIBLE  ORBITING  BODY 
2.1  -  Assumptions 

In  order  to  achieve  a  low  mass  to  area  ratio  many  of 
the  proposed  large  space  structures  have  been  designed  in 
the  form  of  lattice  or  truss  structures.  A  finite  element 
analysis  of  such  an  orbiting  structure  would  require  a  large 
computing  capability  and  may  be  expensive.  A  preliminary 
insight  into  the  dynamics  of  the  system  can  be  obtained  by 
representing  the  structure  as  a  large  thin  plate. 

Early  analyses  of  space  structures  were  based  on 
aluminum  or  aluminum  alloys;  since  then  advances  in 
technology  have  made  composite  graphite  a  feasible 
alternative.  When  comparing  the  two  materials,  graphite 
displays  two  advantages,  flexibility  and  weight,  thus, 
makino  graphite  the  optimum  material. 

The  material  property  v.'.lues  adopted  for  reinforced 
composite  graphite  here  are:  Young's  Modulus,  E,  40  x  10* 
Ib/in^;  Poisson's  ratio,  u,0.3;  and  density,  p,  5,42  x  10^ 
Ib/in^.^^®’  The  structure's  dimensions  are  assumed  as;  width 
and  length,  £,  100  m;  and  thickness,  t,  0.01  m.  It  assumed 
is  to  travel  in  a  circular  orbit  at  an  altitude,  h,  250 
n. miles,  while  maintai  ing  a  constant  angular  velocity,  , 
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0.0011162  rad/sec. 


The  equations  of  motion  are  derived  using  a  fjewton- 
Euler  formulation.  The  principal  assumptions  made  for  this 
development  are: 

(1)  the  mass  is  idealized  as  a  continuous  distribution  of 
mass  points  in  a  plane; 

(2)  the  structural  properties  are  uniformly  distributed; 

(3)  the  material  of  the  body  is  isotropic; 

(4)  the  structure  is  deformed  in  such  a  manner  that  it 
experiences  only  small  strains  (within  the  linear 
range)  ; 

(5)  the  elastic  displacements  are  small  as  compared  with 
the  characteristic  linear  dimensions  of  the  body; 

(6)  the  elastic  deformations  in  the  plane  of  the  plate  are 
much  smaller  in  comparison  to  the  deformations  norrrtal 
to  the  plate; 

(7)  the  first  three  elastic  modes  will  be  considered,  since 
normally  only  a  few  elastic  modes  contribute 
significantly  to  the  vibrational  motion  of  the 
structure ; 

(8)  the  system  is  considered  closed,  i.e.,  no  mass  transfer 
across  the  system  boundaries;  and 

(9)  there  are  no  geometrical  constraints  on  the  motion. 
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2.2  -  Coordinate  Frames 


Initially  the  equations  of  motion  are  derived  for  a 
flexible  orbiting  body  of  arbitrary  shape.  Figure  (1) 
shows  a  flexible  orbiting  body  with  various  symbols  and 
coordinate  frames. 


□  X.Y.  Z.-  Orbit  fixed  frone 


D’  -  center  of  the  earth 
D  -  center  of  nass  of  the  body 


Figure  1.  Transformation  of  Coordinate  Frames 


I?", 

T^:  O'XYZ  is  an  earth  centered  inertial  reference  frame 

(Tj^)  with  O' 2  along  the  earth's  spin  axis  and  O’X  along 
the  ascending  node. 

T^:  Oi^ijij  is  the  local  intrinsic  frame  (t,)  centered  at 

the  center  of  mass  of  the  body  0,  with  Oi,  along  the 
local  vertical  and  Oi^  perpendicular  to  the  plane  20 'O. 

T^:  OXjjYpZ^  is  an  orbit  fixed  reference  frame  (t^)  centered 

at  the  center  of  mass  of  the  body  O,  with  OX^.  along  the 
local  vertical  and  OY^  along  the  orbit  normal  opposite 
to  the  orbit  angular  momentum  vector. 

Tj:  Oxyz  defines  the  principal  axes  of  inertia  (Tj)  of  the 

body  in  the  undeformed  state  (not  shown  in  figure) . 


The  above  reference  frames  are  related  to  each  other  as 
follows; 
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The  various  transformation  matrices  are^"’: 

[sintjcosa)  sirnisincj  cosh  ] 

•  1  f  4  1 

r,  -  jCOSTICOSO)  COSTlSmCi)  -Sir.Tli  '  ' 

[  -sinti)  COSO)  0  J 

T,  represents  the  transformation  from  the  inertial  frame  to 
the  intrinsic  frame. 


10  0 
0  cosx  sinx 
0  -sinx  cosx. 


(5) 


Tj  represents  the  transformation  from  the  intrinsic  frame  to 
the  orbit  fixed  frame. 


cob^cobQ 

-sin^cose 

sin0 


(sin4>cost(r  ■►cos4iein6sinf ) 
(co8  4'cos  <r-sin<>sin6sini|») 
-cosOsiniJt 


{3in4>siniti  -coB^sin6cosf )  j 
(cos4)sint|f +  sin4>sin8co9<i) 
cosScosiJ^  J 


(6) 


Tj  represents  the  transformation  from  the  orbit  frame  to  the 
body  frame  according  to  the  sequence  (1)  \|)-yaw;  (2)  6~pitch; 
and  (3)  0-roll,  respectively.  The  body  angular  velocity 
components  w  ,  u  ,  and  u  are  then  related  to  the  Euler 
angular  rates  0,0,  and  as  follows: 


-  6sin4>  +  i|fcos4>cos0-u)^(sin(|>cosij/  +  cos(t>sin6sinilr)  (7) 

tOy-  0cos4>-Tjrsin4)cos9-w^(cos(J>cosi|/-sin<|)sin0sinttr)  (8) 

o  J- iJrsinO +  (i)  c (cos6sini|f) 


2.3  -  Gravitation 

The  gravitational  potential  at  any  point  can  be 
expressed  in  its  roost  general  forro^®'”-’^^ : 

VCp,!!,^)-  -^^]  +  vay' (  — f  'ft  (ri ,  w)  (10) 

V  P  /  ^  '  P  ' 

where, 

ft(Tl,<*i)  “52  (t))  COS  (m(i> +4>^„)  ]  (11) 

m-0 

represents  the  ro—  associated  Legendre  function  of 
order  s;  and  <t>^^  are  constants  obtained  through  analysis 
of  the  satellite  orbital  motion;  p  represents  the  distance 
from  the  point  to  the  center  of  the  earth;  and  r?  and  u 
represent  the  colatitude  and  the  longitude,  respectively,  of 
the  point. 

For  the  gravitational  force  per  unit  mass  at  the  body's 
center  of  gravity, 


1 


For  a  point  at  a  distance  r  from  G,  neglecting  small 
quantities  of  the  order  jrl/p,  the  gravity  force,  F(x): 


F(x)  r 

Substitution  of  the  matrix  operator,  B*,  leads  to; 


(13) 
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where  and  B^®^  are  defined  as: 


(TjTj)  -’r 
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B**’  - 
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sin'^Ti 


(16) 


where,  (  )'  -  ;  (T)  - 


do) 


Reprojecting  on  the  body  fixed  axes  results  in; 


f{x)  -TjTjFc* 


va‘ 


s) 


(17) 


are  symmetric  matrices,  where  the  following 


and 

definition  of  m”  represents  both  matrices, 

M'-’  -  T.  B"  (T,  (18) 


Affoi 


3cos^4)cos^0 -- 1  -3sin4)cos()>cos^0  3cos(t>cos  0sin0  ] 

-3 sin(t»cos 4)cos^0  3sin-4)Cos^6 -  1  ~3sin<)>ccs0sin0|  ^^8) 
Bcos^'COsSsinO  -3sin(}>cos6sin0  3sin''0-l  j 


-  EE 


B, 


'nk 


Jt-l 


(20) 


where  t^^  is  the  (mn)  element  of  the  (TjT^)’’  matrix  and 
is  the  (mn)  element  of  the  matrix. 


2.4  -  Equations  of  Motion 

The  position  of  a  general  point  with  respect  to  the 
body  fixed  frame,  Tj,  is  given  by 

I  -  q  (21) 

where  r^  represents  the  position  vector  of  the  body  with 
respect  to  O  in  the  undeformed  state,  +  ^ r  and  q 

represents  the  elastic  displacement  of  the  body.  Therefore, 

T-rT  (22a)  and  T-rT  (22b)  (22) 


For  small  amplitude  elastic  displacements,  q  can  be 
represented  as  a  superposition  of  various  modal 
contributions : 


n-1 


(23) 


where  A^(t)  represents  the  n—  modal  amplitude;  and 
represents  the  eigenmodes  of  vibration,  +  j . 

Substitution  of  Eqn.  (23)  into  Eqn.  (17)  forms: 


f-  *y  kJ 


r 

^73-1 

(24) 


The  linear  operator  if'(q]  transforms  small  structural 
displacements,  q,  to  the  structural  forces  acting  on  the 
generic  point  of  the  body.  The  mode  shape,  is 

associated  with  the  natural  frequency,  and  satisfies  the 
following  orthonormality  condition: 

dm  -  5^^  (25) 

where  represents  the  generalized  mass  in  the  nth  mode. 
The  linear  operator  of  becomes; 


dm 


(26) 


Therefore , 


S£[g(r„,  t)  ]  dm-  [-  t  u>lAAt)  ]  dm  (27) 

D'l 

To  satisfy  the  laws  of  conservation,  an  unconstrained  body's 
elastic  inodes  must  be  orthogonal  to  the  rigid  body  modes: 
Conservation  of  Linear  Momentum-Translational  Orthogonality 

j  dm  -  0  (28) 

vol 

Conservation  of  Angular  Momentum-Rotational  Orthogonality 

J  ’"’dm  -  J 

vo2  vol 

If  the  body  is  constrained  against  translational  and 
rotation  of  the  undeformed  center  of  mass,  the  corresponding 
inodes  are  called  "fixed  modes".  For  fixed  modes  Eqns.  (26) 
and  (27)  do  not  hold.  For  fixed  inodes  the  center  of  mass  in 
the  deformed  state  no  longer  coincides  with  the  origin  of 
the  body  frame.  Only  for  the  free  modes  does  Jr  dm  =  0. 

2.4.1  -  Equations  of  Rotational  Motion 

Returning  to  the  development  of  the  local  equations, 
the  equation  of  motion  for  an  elemental  mass,  dm,  whose 
instantaneous  position  vector  from  the  center  of  mass  of  the 
body  is  r,  can  be  written  as: 


-r. 


^2  ] 
0 


dm  -  0 


(29) 
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[cf!  +  f  dn?+ E- adm  (30) 

where  f  represents  the  gravitational  force  per  unit  mass; 

E  represents  the  external  forces  (other  than  gravitational) . 

Rewriting  the  inertial  acceleration  of  a  differential 
mass  as  a  contribution  of  terms  as  seen  by  an  observer  in 
the  rotating  body  fixed  system  of  coordinates,  Tj,  it  can  be 
seen  that  the  general  force  equation  for  a  rotating  body  is: 

adm- dm [?“+ 2cSx2r- tJxr  +  oix  (uSxr)  ]  (31) 

where  3  represents  the  angular  velocity,  o^o^i+u^j+o^R , 

The  equations  of  rotational  motion  of  the  body  are 
obtained  by  talcing  the  moments  of  all  the  external, 
internal,  and  inertial  forces  acting  on  the  body.  After 
equating  Eqns.  (30)  and  (31)  and  talcing  the  moments: 

f rx  [a_ +  ■?*•*■  2 u)x?‘+Tlyxr  +  (3x  (wxiO  ]  pdv-  frx  (-^1^  -t-  /  +  -§-]  dm  (32) 

J  J  am  dm 

The  various  terms  of  Eqn.  (32)  can  be  evaluated  using  vector 

calculus.  Assuming  l4^f|«  1,  only  the  1^  order  terms  in  q 

are  retained.  Through  the  substitution  of  the  values  of  the 

integrals  into  Eqn.  (32)  and  rearrangement  of  terms,  one  can 

obtain  the  following  form  for  the  rotational  equations  of 

motion: 


*T2  **yy 

r  (n)  ,  rr  (n) 


x~  y 

.2  ..2 


(35) 


G)  ,0) ,  )  *  ( w",-G)t)  (/fjf  )  ] 


Qy^^’and  are  obtained  by  the  cyclic  permutation  of  x,  y, 

z  in  the  expression  for  Q/"’,  where  is  defined  as, 


vcl 

The  term,  C,  represents  the  external  momentum  caused  by 
external  forces  other  than  gravitational,  where  e  is  the 
external  force  per  unit  mass. 


C-  Jfxedin 


(37) 


The  terms,  account  for  the  difference  in  position 

between  the  actual  center  of  mass  and  the  center  of  mass  of 


the  undeformed  body.  For  the  case  of  free  modes  0  . 

^  t  dm  (38) 

n-l  J  n-l  J 

The  term,  G^.,  corresponds  to  the  gravitational  torque  on  a 
rigid  body, 

^7-  f  f^xMf^dm  (39) 

voi 

Gr“  ^  ^  z~  ^  y'>  ^22^  *  ^  ^  x~  ^ 

Where  is  an  element  of  the  M  matrix  which  is  defined  as 

(a/p) 

The  terms,  correspond  to  the  gravitational  torque  due 

to  elastic  deformations,  where  G^"^=G  ^"’i+G  ^"^j+G 

X  y  z 

Eg*"’-  f  [f^xMg*qxMrJdm  (41) 

n-l  J  , 

voj 


ry  (n) 

G'x 


-A„[  (W33-M22) 


Ayy 


*XZ 
r(/))  \ 


■M- 
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(42) 


Gy‘"’and  G^^"^  components  are  obtained  by  the  cyclic 
permutation  of  x,  y,  z  in  the  expression  of 
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2.4.2  -  Generic  Mode  Equation 


The  generic  mode  equation  is  obtained  by  taking  the 
modal  components  of  all  internal,  external,  and  inertial 
forces  acting  on  the  body^®-"*’^^ . 


J  $  2o5  xt'+TS'xz  +  Cx  (€)xf~}  ]  pdv 

vol 

f  ^ -[a? [(^ /p  +  f+ ^  pdv 

vol 


(43) 


The  various  terms  appearing  in  Eqn.  (43)  can  be  expanded. 
With  the  substitution  of  the  various  integrals  and 
rearrangement  of  the  terms,  the  following  generic  modal 
equation  results: 


(44) 


m-1 


i7>-l 


in) 


in) 

yz 


^  y  ^  ^xy 


■H, 


In, 

yx 


yO) 


(Ny"  *H: 


T 

in) 

XX 


xz 

in'  , 
zy  ' 

(;3> 

zz 


iH: 


in) 


l^yx  )  ^ 


/  rr  rr 


(45) 


(46) 


) 


The  term,  represents  the  influence  (cr  force)  of  the 

rotational  motion  of  the  body  on  the  n—  elastic  mode  due  to 
inertia . 


^  ^  mn"  f  -tSxq-^  ^  -CO X  {i3x<^  ]  pdv  (47  ) 


(48) 

The  terms,  (p^,  represent  the  influence  of  the  other  elastic 
modes  on  the  n—  elastic  mode  due  to  inertia,  where  is 

defined  as, 

L:r-fK'-K'dK  (45) 

vo} 

The  term,  g^,  is  the  gravitational  force  acting  on  the  n— 
mode  due  to  the  rigid  body  motion. 

f  (50) 

•  -  a  6 

vci 

The  terms,  g^,  represent  the  gravitational  force  acting  on 
the  n—  mode  due  to  the  elastic  motion  on  the  m—  mode. 

tg^-  f  ^‘'^’Mgdw-A^TXLJr'Kii  (51) 

viv  “  » 

The  term,  E^,  is  the  external  force  component  acting  on  the 
n—  mode, 

(52) 

vol 

The  term,  D'^,  represents  the  term  corresponding  to  the 
displacement  of  the  center  of  mass  due  to  the  elastic 
motion.  For  unconstrained  motion  D'  vanishes. 

n 


(p -  2;i  J  ca  J  ^  - 

Ajd),(L;r'-Lir) 

-co^{L;r-Ljr*)  -u)V(Lir-Lr’) 


The  equations  of  translational  motion  of  the  body  in 
orbit  are  obtained  by  integrating  the  sunuTiation  of  Eqns. 

(30)  and  (31). 

J  +  203  xT-tJxr  +  c5x  (coxi^  ]  pdv-  J  1 0  /  p  *  f  *  pdv  (54) 

voj  vol 

Noting  Jrpdv=0,  one  can  obtain  the  translational  equation  of 
motion. 

a^*£^*D-E/w  (55) 

where  is  the  inertial  acceleration  of  the  center  of 
mass;  is  the  intensity  of  the  gravitational  force  at  the 
center  of  mass  of  the  body  (force/mass);  m  is  the  mass  of 
the  body;  E  is  the  resultant  of  the  external  forces;  and 

[;?'+ 2oixg+T!lrxg*t3x  (uxc^  -  S£[q] /p -Wgl  pcfv'  (55) 

VO  I 

For  the  free-free  mode  shapes  -unconstrained  mode  shapes- 
D=0. 

Because  of  the  size  of  the  structure  considered  in  this 
work,  (i.e.,  (100  m) )  one  may  neglect  the  effect  of  elastic 

motion  on  the  orbital  motion  of  the  center  of  mass.  In  this 
thesis,  the  equations  of  orbital  motion  are  not  considered 


for  further  analysis.  It  is  assumed  that  the  orbital 
position  of  the  center  of  mass  can  be  readily  computed  using 
techniques  of  orbital  mechanics. 

In  summary,  the  motion  of  an  arbitrary  flexible  body  in 
orbit  is  described  by  Eqns.  (33),  (44),  and  (54). 

Consideration  of  the  effects  of  gravity  grad^-’nt  -including 
higher  order  harmonics-  has  been  included  in  the  derivation 
of  these  equations.  Equations  (33)  and  (54)  are  vector 
equations  which  describe  the  dynamics  of  the  rigid 
rotational  and  rigid  translational  motions,  respectively. 
Equation  (44)  is  a  scalar  equation  which  describes  the 
elastic  motion  ox  the  body  in  its  n—  elastic  mode.  With  the 
calculation  of  the  natural  frequency  and  modal  shape 
functions,  the  equations  of  motion  of  a  particular  system 


can  be  derived. 


2.5  -  Notion  of  a  Thin  Flexible  Plate  in  Orbit 


One  class  of  structures,  which  has  been  proposed  for 
use  in  many  future  space  applications,  has  the  basic  form  of 
a  thin  flat  plate.  Included  among  the  proposed  applications 
are;  solar  energy  collection,  communications,  and  scientific 
data  based  orbiting  platforms.  A  brief  development  of  a 
mathematical  model  for  the  attitude  motion  and  elastic 
motion  for  a  large,  yet  thin,  flexible,  flat  plate  in  orbit 
is  presented  in  this  section.  The  equations  presented  here 
for  the  plate's  motion  are  the  result  of  simplification  of 
the  general  set  of  equations  that  describe  the  motion  of  an 
arbitrary  flexible  body  in  orbit  (previously  presented  in 
Section  2.4) . 


By  taking  into  consideration  assumption  (9) ,  one  can 
further  simplify  the  development  of  the  equations  of 
rotational  motion.  Assumption  (9)  leads  to  where 

for  all  a  and  0,  similarly,  D^=Dj^j=0 ,  With  a 

rearrangement  of  terms  and  the  separation  of  the  rotational 
equation  of  motion  into  its  vector  components,  Eqn.  (33)  can 
be  developed  to  yield  the  following  set  of  rotational 
equations  of  motion  for  the  elastic  plate  in  orbit. 


ly^  y  {  I ^  I y)  ^  -W  y  ~  Cy  *  G j.y 

J  (  ly~  ly)  U)  y(t>  ^  ^  ^ ^ 


(57a) 

(57b) 

(57c) 


(57) 
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The  general  equation  for  elastic  motion  is  listed  below. 

m-l  m-1 

This  paper  will  present  three  different  nominal 
orientations  of  the  platform  in  orbit;  attitude  and  shape 
control  will  be  achieved  for  each.  They  are: 

Case  (1)  the  platform  following  the  local  horizontal  with 
its  larger  surface  normal  to  the  local  vertical, 
see  Figure  2.1; 

Case  (2)  the  platform  following  the  local  vertical  with  its 
larger  surface  perpendicular  to  the  plane  of  the 
orbit,  see  Figure  2.2;  and 

Case  (3)  the  platform  following  the  local  vertical  with  its 
larger  surface  perpendicular  to  the  orbit  normal, 
see  Figure  2.3. 

It  can  be  shown  that  for  gravitational  stability  the 
plate's  axis  of  minimum  moment  of  inertia  should  be 
nominally  oriented  along  the  local  vertical.  However,  in 
many  applications  it  is  required  that  the  major  surface  of 
the  plate  be  pointed  towards  the  earth.  Therefore,  fo'  the 
plate  orientation  of  Case  (1)  a  complete  development  of  the 
equations  will  be  presented.  For  the  two  other  orientations, 
Case  (2)  and  Case  (3),  the  final  form  of  the  equations  of 
motion  will  be  shown. 


Figure  2.1.  Case  (1)  •>  Platform  along  local  horizontal 


Figure  2.3.  Case  (3)  Platform  following  local  vertical  with 
major  surface  in  the  orbit  plane 
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2.5.1  -  Plate  Normal  Nominally  Along  the  Local  Vertical 

Case  (1) 

From  the  equations  of  motion,  Eqns.  {57a, b,c),  and  (58) 
can  be  written, 


Yaw: 

(59a) 

J 

-^x 

Pi tch: 

(Jx--fz)l 

(59b) 

Roll'. 

(ly-I,)  ] 

I. 

(i>  Ct>  * 
y'"  X 

I. 

(59c) 

(59) 


Generic 

Mode 


m-i 


m-1 


For  Case  (1)  the  plate  is  oriented  with  its  normal 
following  the  local  vertical,  its  axis  is  aligned  with  the 
outward  direction  of  the  local  vertical  (see  Fig.  2.1). 

Under  the  nominal  motion  the  Euler  angles,  i)i,  <t>,  and  6,  are 
defined  according  to  the  sequence  given  after  Eqn.  (6).  From 
the  previously  shown  transformations,  the  Euler  angular 
rates  are  related  to  the  body  rates  as  derived  in  Eqns.  (7) , 
(8) ,  and  (9) .  In  order  to  examine  the  stability  of  the 
system  Eqns.  (59)  can  be  linearized  by  assuming  small 
amplitude  pitch,  roll,  and  yaw  displacements  and  also  small 
values  of  their  respective  time  derivatives.  With  this  in 
mind  the  angular  rates  become: 


cJjj-4f-Ci)^<J)  —  (61) 

ojj.  -  6-  Ci)^  —  (0^,-6  (62) 

G)^-4)+tJ^4'  “•  6>jp-(i>+(i)^ijf  (63) 

As  a  result  of  the  previous  assumptions,  (2),  (3),  (6),  (7) 

and  (9) ,  further  simplifications  of  the  equations  of  motion 
are  in  order.  Assumption  (6)  leads  to  (r^)  , 

where  n  is  unit  normal  vector  to  the  plate.  As  previously 
stated,  by  assumption  (9),  and  and 

^n=D(„)=0.  In  addition,  where, 

S^-Kroneckez  delta  (64) 

Hence, 

g‘"''-0;  tp„-0;  g^-0  (65) 

<p  im,  -  -  (  “  y  -  W  -  (  6)  y-  6)^,)  6  (  66  ) 

G..  -  (  I.  -  ly)  ;  G-y  '  U,  -  ^31  ;  ^1.  (  68  > 

Assumptions  of  a  spherical  symmetric  gravitational 
field  and  circular  orbit  and  neglecting  the  higher  order 
terms,  ,  of  the  M  matrix  results  in  the  following 

simplification  of  M23,  and 


n  /  *>  2 

va  “  /  p  -  o)  c 


-3u)  cSin(}>cos4)cos^6 


'  <  O  ) 


Af-,  -  -3co“,-sin4'cros0  sin6 

(70) 

-  3a)^,-cos4)Cos6sin6 

(71) 

with  the  use  of  the  above  simplifications  -for  this 
case-  the  yaw,  roll,  pitch,  and  the  generic  mode  equations, 
respectively,  become. 


[3sin<))cos6sin6] 


[$<{;- fa)  ^(<{;-t-64))  +U)e4>] 

^  z 

Sl  -  f3sin(}>cos4>cos^6] 

•^z 


(73) 


©  +  u)  ^(i)>\ii-<i>4>)  -  o'pijjt})]  - 


^  SIjljLiI. [3cos(}>cos6sin6] 


(74) 


(75) 


where  is  defined  as 


Af,!  -  Wc  [3cos^4>cos^8-l] 


(75 
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The  equations  of  rotational  motion  can  be  simplified  further 
with  the  elimination  of  higher  order  nonlinear  terms 
involving  the  products  of  two  angles  or  angular  rates  or 
products  of  an  angle  with  an  angular  rate,  and  assuming  the 
rigid  body  angular  oscillations  are  small,  i.e., 

The  resulting  equations  are  presented  below: 


- 

1 - 1 _ * 

w  C 

I, 

1,-1, 


—iL - 1  14  4). 
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I, 


(77) 


(78) 
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(3  <*>^,,6)  -Cy/ly 


(79) 


At  this  point,  further  development  of  the  generic  mode 
equation  may  proceed  with  the  substitution  of  the 
appropriate  terms.  This  leads  to 


[Wn-  (6'^- 20)  ,-6  +  +  *2(i)  +  U) 


U)  ^  (  3  005-4)005^6-1)  ] 


(80) 


Assuming  the  Euler  angle  displacements  and  time  derivatives 
and  the  transverse  elastic  deformation  amplitudes  and  time 
derivatives  are  small,  with  the  removal  of  nonlinear  terms 
Eqn.  (80)  simplifies  further  to, 


u> 


(81) 


Based  on  the  preceding  assumptions  the  following  linearized 
equations  of  motion  are  obtained  in  a  dimensionless  form; 


T  -  U> 

:  ;  4>^-d<|>/a’T  ;  etc. 

Yaw. 

+ 

(82) 

Roll: 

4)"  +  ,4>  +  ( 1  -  Q  -  (C.*  TJ  !  1^,0)^ 

(83) 

Pi  tch : 

^'-■h^y^-iCy^Ty)  llyO\ 

(84) 

Elastic  . 
modes 

(85) 

where  is  defined  as, 

Qj,-  (Ij-Iy) /Ij, ;  Qy-  fly  i  a  ily- I^) /I^  ; 

T  represents  the  torque  produced  by  the  actuators,  and 
represents  the  generic  force  on  the  r—  mode. 

Assumption  (7)  should  be  noted  here,  i.e.,  only  the 
first  three  modes  are  considered.  The  following  observations 
can  be  made  from  Eqns.  (77-80)  regarding  the  motion  of  a 
flat  plate: 

a)  Uncontrolled  rigid  body  motion  is  independent  of  the 
uncontrolled  elastic  motion  of  the  plate. 

b)  The  uncontrolled  elastic  modes  are  coupled  to  the  rigid 
body  motion  through  both  inertia  and  gravity. 

c)  To  first  order,  gravity  can  not  excite  the  elastic 


modes . 


d)  There  is  no  intercoupling  between  the  uncontrolled 


elastic  modes 

either  through  inertia  or 

gravity . 

For  Case  ( 1) , 

I  =21  =21  .  the  rotational 

X  y  z 

equations  of 

motion  become. 

(86) 

4>"-44)*2i(r'-  (C^+  rj 

(87) 

(Cy^Ty)  /I 

(88) 

Since  there  are  no  inputs  from  the  elastic  terms  in 
Eqns.  (82-84)  the  uncontrolled  rigid  body  rotawional  inodes 
are  not  influenced  by  the  elastic  motion  of  the  plate. 
However,  the  elastic  modes  are  coupled  to  the  rigid  body 
rotational  modes  through  higher  order  nonlinear  terms  as 
shown  in  Eqns,  (77-80) .  Further,  for  this  orientation  of  the 
plate  ny>0,  the  pitch  motion  is  unstable  in  the  absence  of 
external  restoring  torques,  (note  the  term,  in 

Eqn.(84)).  For  the  class  of  large  space  structures 
considered,  passive  control  has  been  analyzed,  (see  Kumar, 
ref.  8),  but  for  this  thesis  active  control  will  be  applied. 
Through  the  application  of  the  LQR  problem  with  discretized 
input  data  and  point  actuators  located  on  the  plate's 
surface,  active  control  will  be  implemented. 
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2.5.2  -  Plate  Normal  Nominally  Along  the  Local  Horizontal 

Case  ( 2 ) 

Assumption  (6)  in  this  case  leads  to  (r^) (h)  = 

♦  Similar  to  Case  (1)  ,  one  can  show  that  for  this 

orientation: 

^<">-0;  (p„-0;  g„-0; 

‘P/nn--  and, 

in  Eqns.  (33)  and  (44).  Thus,  for  Case  (2)  the  equations  of 
rotational  motion  are  exactly  the  same  as  Eqns.  (82-84), 
with  the  exception  that  the  relationship  between  the 
principal  moments  of  inertia  and  I,,  differs, 

Ij®=2Ij^=2Iy.  Hence,  the  Yaw,  Roll,  and  Pitch  equations  of 
motion  become,  respectively. 


(89) 

4>"-V  -  /I 

(90) 

e''-3e-  iCy^Ty) 

(91) 

The  equations  of  elastic  motion  are  based  on  Eqn.  (60) , 
together  with  the  above  assumptions  there  results. 


(92) 


where  Mjj  =  ( 3sin^0-l )  .  Removing  angular  crossterms  a*jJ 

restricting  the  higher  order  terms  of  the  angles  yields, 

[Q^„-O)^,-(<j)"-3  0')  +2u)^e-  (93) 

For  small  amplitude  angles  and  elastic  motions,  the  non- 
dimensionalized  the  generic  mode  equation  can  be  written. 


(94) 


Notice  for  small  pitch  amplitude,  Eqn.  (93)  can  be 
reduced  to  a  Mathieu  type  equation.  The  roll  and  yaw 
motions  are  gyroscopically  coupled  to  each  other  and  they 
are  not  influenced  by  either  the  pitch  or  elastic  modes 
within  the  linear  range. 


2.5.3  -  Plate  Normal  Nominally  Along  the  Orbit  Normal 

Case  (3) 

In  this  last  case,  assumption  (6)  leads  to 
(tj,)  (h)  •  Following  the  steps  of  the  derivation 

of  the  two  previous  Cases,  (1)  and  (2),  one  will  arrive  with 
the  same  equations  of  rotational  motion  as  Eqns.  (82),  (83), 
and  (84),  with  the  difference  that  the  principal  moments  of 
inertia  I^,  and  1^,  -for  Case  (3)-  have  the  relationship 
Iy=2Ij,=2I^.  Hence,  the  rotational  equations  of  motion  become; 


4  1 


(95) 

<t)^'*44)-  !C,+  r.)  //,wV 

(96) 

(97) 

The  equations  of  elastic  motion  are  now, 

(98) 

where  ~  (^sin^(pcos^B-l)  .  Cancellation  of  the  angular 
crossterms  and  restriction  of  the  motion  to  only  small 
amplitudinal  pitch  motion  produces. 

Removing  higher  order  and  nonlinear  terms,  plus  non¬ 
dimens  ionalizing  the  amplitude  leads  to, 

(100) 

In  this  final  case,  roll,  pitch,  yaw,  and  the  generic 
modes  are  decoupled  from  each  other.  Notice  without  ary 
external  influences  the  generic  modes,  yaw,  and  pitch  motion 
exhibit  simple  harmonic  motion  while  the  pitch  rate 
increases  linearly  with  time  for  a  given  initial  pitch  rate. 
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-  DETERMINATION  OF  MODAL  FREQUENCIES  AND  KODK 
SHAPES  FOR  A  THIN  FLEXIBLE  FREE-FREE  PLA T E 

The  ability  to  accurately  determine  the  frequencies  and 
mode  shapes  is  essential  for  the  analysis  and  control  of 
large  space  structures  in  orbit.  Four  different  methods  have 
been  analyzed : 

1)  The  approximate  frequencies  and  mode  shapes  of  a 
rectangular  plate  are  obtained  from  a  Rayleigh-Ri tz 
formulation  by  Warburton ; 

2)  The  analytical  results  for  a  square  plate  are  derived 
from  the  Rayle igh-Ritz  method  by  Lenke;'^^^ 

3)  The  frequencies  and  mode  shapes  are  computed  using  the 
finite  element  program,  STRUDL,  developed  at 
Massachusett-s  Institute  of  Technology  r  and 

4)  The  rigid  body  modes  and  the  mode  deflection  are 
calculated  with  an  updated  version  of  STRUDL,  GTSTRUDL, 
developed  at  Georgia  Institute  of  Technology . 

A  tabulation  of  the  numerical  results  for  a  free-free 
thin  aluminum  square  plate  with  the  length  and  width  of 
100.0  m,  thickness  of  0.01  m.  Young's  Modulus  of  0.7441  x 
lo’*’  kg/m^,  Poisson's  ratio  of  0.33,  and  mass  density  of 
2768.0  kg/m^'^’’^ 
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3.1  -  Formulation  bv  Warburton 


An  application  of  the  Rayleigh  Method  allows  the 
derivation  of  an  approximate  frequency  formula .  The 

plate  vibrational  equation  in  the  cartesian  coordinate 
system  (x,y)  is  used  as  a  basic  ecfuation,  with  the  length 
and  width  of  the  plate  taken  along  the  x  and  y  directions, 
respectively,  and  is  given  as. 


,  I2p(l-v") 

dx*  dx^dy^  dy*  Egh^ 

where  p,  v,  and  E  are  the  weight  density,  Poisson's  ratio, 
and  Young's  modulus  of  the  plate  material,  respectively;  h 
is  the  plate  thickness;  and  g  is  the  acceleration  due  to 
gravity.  The  displacement,  w,  is  a  sinusoidal  function, 
which  at  any  point  (x,y)  at  any  time  t,  is  given  by, 


w- iVsinti)  t- a6  (x)  4)  ( y)  sinu) t  (102) 

where  0(x)  and  4>(y)  are  beam  functions  orthogonal  to  each 
other  and  arc  used  to  approximate  the  plate's  behavior;  and 
<i>  is  the  vibrational  circular  frequency . The  appropriate 
free-free  beam  functions,  6(x),  are  defined  as. 


0(x)-l  for  m-O  (103a) 

e(x)-l-~  for  w-1  (103ib) 

a 

B  (x)  kcoshy^-^-—^  for /n- 2 , 4 , 6  , (103c) 

6(x)  -siny''/— -  — ]  +  k'sinhy'|— j  for  /n-3,5,7,... 


(103) 


(103d) 


where , 
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sir.ji- 
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1  1 

can  —  Y  tanh  —  y  -  0 

2  ^  2 

(104a) 

sinli  -i- Y 

2 

sin --y' 

2 

and 

tan  —  Y^ "  tanh  -i  y'^  "  0 

2  2 

(104Jb) 

sinh  —  Y'' 

2 

(104) 


The  corresponding  expressions  for  0(y)  can  be  obtained  by 
substituting  y,  b,  e,  and  c  for  x,  a,  y,  and  k, 
respectively.  After  calculating  6(x)  and  0(y),  the  frequency 
expressions  for  a  rectangular  plate  are  derived 

,  ;  _  paM2Tif)^12  (1  (105) 

n"*  Egh^ 

+  +  2~  [vH^Hy-^  (1-v)  J^y]  (106) 

where  1  is  a  non-dimensional  factor,  proportional  to  the 
frequency;  a  and  b  are  the  length  and  width  of  the  plate; 
and  G,,  H,,  J,,  G^,  and  are  functions  associated  with 

the  number  of  nodal  lines  (m  and  n,  parallel  to  x  and  y)  and 
the  boundary  conditions.  From  Eqn.  (105),  the  frequency  is 
obtained  as: 


IhTC 

Eg 

48p  (l-v2) 

1/2 


(107) 


This  frequency  expression  is  valid  for  thin  rectangular 
plates;  for  square  plates  it  must  be  modified.  For  square 
plates,  the  combinations  of  (m,n)±{n,m)  describe  the  types 
of  existing  modes  (see  Fig.  3.1-3.10);  for  these  cases  k  in 
Eqn.  (107)  must  be  modified. 

Modes  fm.01±(0.m^  .  for  m  =  2.4.6... 


A.2-  (m-Jl)*±2v  (/n-.l)2. 


Modes  (m.l)±(l.m).  for  m  =  3.5.7, 


(108) 


12.  (1-v)  (;n--i)2 

2  2 


±2v(m-4)"-^ 
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1- 


{M-  )  It 

2 

±2 (1-v) 
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192 


(109) 


For  any  mode  of  vibration  the  nodal  pattern  is  defined 
by  m  and  n  -the  number  of  nodal  lines  in  the  x  and  y 
directions-  respectively.  The  mode  shapes  are  obtained  by 
using  the  corresponding  modal  frequencies  in  the  beam 
functions  and  then  evaluating  the  product,  6(x)-<j(y), 
numerically. 

Using  Warburton's  results,  Eqns.  (108)  and  (109)  and 
the  expressions  for  0(x)  and  <p(y)  ,  frequencies  and  mode 
shapes  are  calculated  for  different  combinations  of  the 
number  of  nodal  lines,  m  and  n,  starting  with  combinations 
of  m=0  and  n=l,  through  m=3  and  n=3 .  The  first  three 
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combinations  of  nodal  line  numbers,  (0,0),  (1,0),  and  (0,1), 

represent  rigid  body  motion.  The  combination  of  m=l  and  n=l 
produce  the  first  fundamental  flexural  frequency.  The 
corresponding  mode  shape  for  the  plate  is  obtained  by 
multiplying  the  beam  functions  6(x)  and  <^{y),  for  the  mode 
characterized  by  m=l  and  n=l  (see  Fig.  3,1).  Since  the  plate 
is  approximated  by  sets  of  orthogonal  beam  functions  in  the 
X  and  y  direction,  the  nodal  pattern  is  also  obtained  by 
plotting  the  nodal  points  of  these  beams  for  their  first 
several  modes.  The  next  two  higher  frequencies  are  obtained 
by  combinations  of  m=0  and  n=2,  but  the  nodal  patterns  can 
not  be  visualized.  This  is  because  the  frequencies  are  of  a 
special  type  resulting  from  a  combination  of  the  (2,0)  and 
(0,2)  plate  modes.  When  the  mode  corresponding  to  (2,0) 

(Fig.  3.2)  is  superimposed  with  -(0,2),  (Fig.  3.3),  the  mode 
shape  depicted  in  Figure  3.4  results.  By  imposing  the  (2,0) 
and  (0,2)  combinations,  the  third  mode  shape  (Fig. 3. 5)  is 
obtained.  The  combinations  of  nodal  patterns  m=l  and  n=2, 
give  identical  frequencies  for  the  fourth  and  fifth  mode  and 
the  corresponding  shape  (Fig.  3.6)  is  as  expected.  The 
following  higher  two  frequencies  are  also  identical  and  are 
the  results  of  combinations  from  the  (3,0)  and  (0,3)  nodal 
pattern  lines  (Fig.  3.7).  The  eighth  frequency  is  obtained 
from  ro=2  and  n=2  and  the  mode  shape  obtained  is  shown  in 
Figure  3.8.  The  higher  frequencies  are  obtained  in  a  similar 
manner.  The  ninth  and  tenth  modes  are  the  special  type 


resulting  from  combinations  of  (3,1)  and  (1,3),  (see  Figs 
3.9  and  3.10)  and  the  next  higher  frequencies  are 
combinations  of  (3,2)  (2,3),  and  (3,3)  nodal  lines, 
respectively.  The  calculated  modal  frequencies  and  nodal 
patterns  are  shown  in  Table 


Figure  3.1.  First  Mode  (1,1)  Figure  3.2.  (2,0)  Mode 


Figure  3.3.  (o,2)  Mode  Figure  3.4.  Second  Mode 


(2,0)-(0,2) 
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Figure  3.5.  Third  Mode  Figure  3.6.  Fourth  Mode 

(2,0)+(0,2>  (2,1) 


Figure  3.7.  sixth  Mode  (3,0)  Figure  3.8.  Eighth  Mode 


(2,2) 


Figure  3.9.  Ninth  Mode 
(1,3)-(3,1) 


Figure  3.10.  Tenth  Mode 
(l,3)+(3,l) 
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3.2  -  Formulation  bv  Lemke 

The  analytical  results  for  a  square  plate  are  derived 
from  a  method  of  Lemke,  which  also  uses  the  Rayleigh-Ritz 
method.  The  displacement  function  w,  can  be  expressed  as 
a  sinusoidal  function  (at  any  point  (x,y)  and  time  t) , 

w-Wsinu)t  (110) 

where  w  is  the  circular  frecpaency  (o=2?rf)  .  The  amplitude  W, 
can  be  expressed  by, 

h'(x,P)  -EA^X^ixl  yjy)  (111) 

where  X^(x)  and  Y^(y)  are  the  free-free  beam  functions 
expressed  in  terms  of  a  normalized  (i.e.,  x=x/a,  where  a=£) 
xy  coordinate  system  having  the  origin  at  the  plate  center; 

,  coshJc-Cos Jc_x+ cosJc„coshJc„x  , 

X„(jc) - - .  -  - - ^  (m  even)  ;  (ii2j 

^cosh^k„  +  cos^  Jc„ 
sinhk_sinic_x+sinjc-sinhk_x  , 

- - £L- - (modd);  (113) 

y/sinh^ic„-sin2ir„ 

Y^(y}  is  obtained  by  replacing  x  by  y  and  m  by  n,  in  Egns. 
(112)  and  (113).  The  values  k„  are  the  roots  of  the 
characteristic  equations,  that  are, 

tanJc„+ tanhk^- 0  [m  even)  (114) 
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tanic„- tanhJc^- 0  {m  odd)  (115) 

These  equations  result  from  the  spatial  boundary  conditions 
For  a  rectangular  plate  the  potential  energy  due  to  bending 
is  given  by. 


Eh^ 

12 


(1-v) 


\2 

idxdy 


(116) 


where  u,  E,  and  h,  are  the  Poisson's  ratio.  Young's  modulus 
and  the  thic)cness  of  the  plate  material,  respectively.  The 
Jcinetic  energy  is  represented  as. 


(117) 


where  p  is  the  density  of  the  plate  material,  and  g  is  the 
acceleration  due  to  gravity.  Thus, 


2_ 

2 


(118) 


^  0  0 

where  is  the  maximum  potential  energy  due  to  bending, 
and  T^^  is  the  maximum  )cinetic  energy  due  to  bending. 

By  setting  =  T^^,  it  can  be  shown, 

- - - - 

^ff«‘dxcly 

'00 


(120) 


The  coefficients  A^,  in  Eqn.  (Ill)  are  determined  to 
produce  in  Eqn,  (120)  -a  minimum,  Lemke  obtained  the 
coefficients  A^,  by  taking  six  or  more  terms  in  the  series, 
Eqn,  (111),  and  using  four  different  values  of  Poisson's 
ratio.  Expressions  for  six  mode  shapes  and  frequencies  along 
with  the  coefficients  A„,  are  tabulated  in  reference  24.  As 
an  example  the  expression  incorporating  15  terms  for  the 
first  mode  is  given  here: 


W{x,y^  +  0 . 0325  -0 . 005X3^3 -_0_.  00257 

+  0.00121  [X^Y^*X^Y^)  -0.000365X575+  .  .  . 


and  (ji  - 


13 . 086 


Eh- 


a2  \  12p(l-v2) 


for  V  -  0 . 34  3 


Frequencies  and  mode  shapes  for  the  first  six  modes  were 
obtained  using  the  Lemke  method  (see  Table  1).^”^ 
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3.3  -  STRUDL,  Structural  Design  Language 

A  computer  program  implemented  was  E .  IJDL  (Structural 
Design  Language It  was  developed  by  the  Civil 
Engineering  Department  of  the  Massachusetts  Institute  of 
Technology  and  installed  on  Howard  University's  mainframe 
computer  (IBM  3090) .  STRUDL  has  the  capability  to  apply  the 
finite  element  method  to  determine  the  mode  shapes  and  the 
natural  frequencies  of  vibration.  Necessary  computer  input 
by  the  user  is  the  specification  of  the  structure's  material 
properties,  dimensions,  and  the  types  and  total  number  of 
finite  elements  representing  the  structure. 

For  a  rectangular  plate,  the  finite  elements  can  be 
specified  as  rectangular  elements;  the  number  of  elements 
into  which  the  plate  is  divided  depends  on  the  accuracy 
required.  In  the  finite  element  method,  the  discretization 
of  the  number  of  degrees  of  freedom  requires  the 
introduction  of  simplifying  assumptions  in  the  element 
formulation,  which  represents  an  important  source  of  error 
in  the  results.  As  a  consequence,  finite  element  results  are 
dependent  upon  the  number  of  elements  used  in  the  model. 
Thus,  even  when  consistent  element  formulat  .ons  insure 
convergence  of  the  results  as  the  number  of  elements  is 
increased,  finite  element  models  cannot  be  arbitrarily 
designed.  Generally,  the  accuracy  of  the  results  will 
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improve  when  the  structure  is  modeled  with  a  higher  number 
of  elements.  However,  computational  errors  due  to  truncation 
and  round-off  errors  may  predominate  as  the  order  of  the 
elements  increases  beyond  a  limit.  Further,  hardware 
calculation  limitations  can  restrict  the  number  of  elements 
into  which  a  structure  is  divided,  thus,  limiting  the 
ability  to  obtain  better  results. 

For  the  modes  specified,  STRUDL  calculates  the 
deflections  at  each  element's  corners;  from  the  deflections 
the  modal  shapes  can  be  determined.  For  each  mode,  a  set  of 
frequencies  is  also  generated.  Results  obtained  from  STRUDL 
are  available  in  Table 
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3.4  -  GTSTRUDL,  Georgia  Tech  Structural  Design  Language 

GTSTRUDL  is  an  updated  version  of  STRUDL  developed  at 
Georgia  Institute  of  Technology .  Similar  to  STRUDL, 
GTSTRUDL,  uses  the  finite  element  method  to  determine  the 
modal  shapes  and  the  vibrational  frequencies.  The 
structure's  material  properties  and  its  dimensions  are 
necessary  initial  inputs.  In  addition,  the  number  of 
elements  into  which  the  structure  is  divided  and  the  type  of 
element  must  be  specified  by  the  user.  For  this  thesis's 
structure,  a  thin  flat  square  plate,  the  'BPR'  element  - 
bending  plate  rectangular  with  square  dimensions-  was 
chosen.  This  element  only  considers  the  deflections 

normal  to  the  major  surface  and  the  in-plane  rotations  at 
the  nodes,  resulting  in  three  degrees  of  freedom  at  each  of 
the  four  nodes  (see  Figure  4). 


In  addition,  the  user  must  specify  the  mass  type: 
consistent  or  lumped.  When  specified,  the  computer  will 
internally  generate  the  mass  matrix.  A  lumped  mass  approach 
will  produce  a  composition  consisting  of  a  diagonal  mass 


matrix,  whose  elements  correspond  to  each  grid  point's 
(node)  mass.  By  evenly  dividing  the  structure's  mass  by  the 
number  of  elements,  and  then  dividing  that  elemental  mass 
into  four,  each  node  is  assigned  some  mass.  If  a  node  is 
shared  by  additional  elements,  its  mass  is  the  total  sum  of 
each  partitioned  mass.  The  total  mass  at  each  node  becomes: 

a)  the  mass  of  the  nodes  at  the  corners,  m^ , 

m,  =  Mj  H-  (#  of  elements  x  4)  ; 

b)  the  mass  of  the  nodes  along  the  sides,  m^, 

m^,  =  (?  of  elements  x  2)  ;  and 

c)  the  mass  of  the  nodes  insioe  the  boundaries,  m^^^, 

{»  of  elements)  . 

3.4.3  -  The  Stiffness  Matrix.  TKl 


The  user  may  input  the  system's  stiffness  matrix  or  it 
may  be  generated  by  the  computer,  the  latter  option  was 
chosen  for  this  thssi.s.  To  begin  the  tormulation  of  the 
stiffness  matrix,  the  stress  (a)  and  strain  (c)  relationship 
becomes  useful. 

' 

dx^ 

^  A  (125) 
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The  generalized  strain  displacement  matrix,  [B],  gives 


the  strains  at  any  point  within  the  element  due  to  unit 
values  of  nodal  displacements . 

{c}  -  (B]  {a}  -  [HI  [«] -MaJ  where  [B]-[H][A]-^  U26) 


{e}  -  [H]  {  a  }  -  [B]  [^J  {a} 


(127) 


From  Hooke’s  Law, 


n 


2 (1-v) 


[  dxdyj  "1  dx~ 


f dxdy-  {a}  ~{g] 


(128) 


Taking  the  variation,  where  (5'''^n=0,  results  in 


//H, 


d^W 

Bx^ 


r  dxdy) 


dxdy-  {6a}  ’’{g)  - 0 


(129) 


After  substituting  in  the  variation  of  the  stress-strain, 
(also  known  as  internal  virtual  work)  the  following  is 
obtained, 


JJ{6e}  ^{oldxdy- {6a}  ’■{g} -0  (130) 

K 

Remembering  the  relationship  ( 6 e  } {  6a )  ^  [  B]  ^  and  substituting 
into  Eqn.  (130),  results  in: 

|j’{6a}  ^[B]  ^[D]  [t]dxdy- {6a}  ^{g}  -0  (131) 


I 
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{6a}  [5]  ^[I>]  IB]  {a}  dxdy~  {q]  ]  -  {6a}  f/fj  {a}  -  {g}  )  -  o  (i32) 

R 

where  {6a)^  is  arbitrary,  and  the  stiffness  matrix,  [K],  is 
defined  as, 

[K]  -  j  j  [3]  ^[D]  [B]  dxdy  ..  [ifl{a}-{g}  (133) 

where  the  nodal  force  vector,  (q),  is  defined  as, 

{g}-/Jp(x,y)  IW}  ‘dxdy  ;  (134) 

and  the  interpolation  matrix,  [N],  is  represented  by  the 
relationship , 

[W3  -  [J’J 

where  the  matrix,  [J] ,  is  defined  as, 

x^.  y^,  xj,  Xjy^,  yj,  x] ,  xjy^,  x^y],  y] .  x^y^.  x^yj]  (136) 

For  a  complete  listing  of  the  stiffness  matrix  individual 
elements  refer  to  Appendix  A. 

3.4.4  -  Development  of  The  Basic  Dynamic  Matrix  Ecmation 

D'Aleiataert 's  principal  incorporated  with  the  principal 
of  virtual  wor)c  (with  inertia  as  a  body  force  distributed) 
can  be  expressed  as^^” : 

-///p  { 6u}  '{ u}  dv+  ///  {6  u}  ~ { b]  dv*  i 6u]  ^{T]  dA~ !!!  [6e  ]  "[  o]  oV-  6  L-  ( 137 ) 

V  V  s  y 
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where  (T)  and  (U)  are  the  kinetic  and  potential  energy,  and 
(u)  represents  the  displacement  field  {u}={u,  v,  w!%  which 
can  be  expressed  in  terms  of  the  nodal  displacements,  (a), 
and  the  interpolation  function,  [N] ,  as  follows, 

{u}  -  im  {a}  (138) 

Substitution  of  [N]  and  {a}  leads  to: 

{u}  -  [J]  {o}  -  [J]  {a)  (139) 

Using  the  interpolation  functions  expressed  in  terms  of 
the  spatial  coordinates  (in  the  static  form)  and  considering 
nodal  displacements,  (a),  as  functions  of  time,  discretizing 
the  domain,  and  making  the  various  substitutions  of  Eqns. 
(123),  (133),  (138),  and  (139)  into  Eqn. (135) ,  forms, 

-fflp{t>a]  {S]dv-^!!!{ba]  ^(JV)  ^{b]dv*H{ba)  '^{r}dA 

V  V  B  (140) 

-///(6a]  ^[iV]  ’’[J]  [1^  {a}dv-///{6a}  ^[B]  ^[D]  [fi]  {a}dv 

V  V 

where , 

Mass  Matrix  ;  [M\  - Hfp  [N]  ‘^[M\  dv  (141) 

V 

Damping  Matrix  ;  [C]  -/// [JV] [W]  dv  (142) 

V 

Stiffness  Matrix:  [K]  -///[B]  '^[D]  [Bj  dv 

V 


(143) 
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Force  Vector  :  {q}  -  !!f  ‘^{b]  dv-^  H  [N]  ‘^{T]  dA  (144) 

V'  s 

Equation  (144)  now  becomes, 

{&a}  (-  m  {S}  4  {g}  -  Id  {d] )  -  {6a]  ^{K]  {a}  <14S) 

After  rearrangement  of  terms,  the  basic  dynamic  matrix 
equation  is  formed; 

[W]  {^}  +  Id  {^}  +  (a)  -  (g)  (146) 


3.4.5  -  Natural  Frecmencv  Formulation  for  Free  Vibrational 
Motion 

Consider  the  free  vibrational  -  undamped  harmonic  - 
motion  system,  (i.e.,  no  damping  and  no  external  forces)  and 
notice  the  damping  matrix  and  the  force  vector  are  null. 
Thus, 

[Af]  ^  [An  {al  -0  (147) 

The  global  equation  for  free  vibration  is 

{a}-{^e-‘“'^  (148) 

where  (a)  is  the  modal  vector  of  constant  amplitude;  and  u 
is  the  natural  frequency.  Substituting  the  global  equation 
into  the  dynamic  equation  results  in. 


( [d  -uMw] )  {^}  -0 


(149) 
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Equation  (149)  is  in  the  form  of  a  basic  eigenvalue- 
eigenvector  problem.  For  a  non-trivial  solution,  the 
determinant  of  the  coefficients  of  {a}  must  be  equal  to 
zero. 

-u)MM]i-o  (150) 

Once  the  characteristic  equation  is  found,  the  roots 
may  be  obtained  by  solving  for  the  eigenvalues  of  the 

n—  mode,  the  natural  frequency,  For  each  natural 

frequency  generated,  an  eigenvector  is  calculated;  these 
provide  the  natural  mode  shapes.  To  ensure  that  the 
solutions  are  linearly  independent,  each  eigenvector  must 
satisfy  the  orthogonality  condition: 

JtAfl  {<?}  j-0  for  i*j  (151) 

Rewriting  Eqn.  (147)  leads  to: 

[in  {5}  -  [if]  {s}  {s)  ^  (152) 

After  multiplication  by  {a}^-  and  {a)^  and  subtraction, 
there  results. 


J[in  {51 I[in  {s]  j[wi  [s]  l[ifi  {51  j  (153) 

It  is  remembered  that  [K]  and  [M]  are  symmetric,  therefore, 

(u"^-o)',)  (J-)  J[in  {5}  -0  (154) 


since , 


6  2 

0  then  {a}  J[Afl  {a}  j  -  0  (155) 

When  orthogonality  is  satisfied  the  roots  are  said  to 
be  distinct.  An  orthogonality  chec)c  was  performed  by 
GTSTRUDL  for  each  eigenvector  solution.  It  basically 
determines  whether  the  normalized  spectral  matrix 
(diagonalized  eigenvalues)  adequately  represents  the 
original  spectral  matrix. 

It  frequently  happens  in  complex  systems  that  the 
frequencies  are  "closely  spaced"  frequencies,  that  is,  cases 
in  which  and  differ  by  only  1%  or  so.  It  occasionally 
happens  that  a  system  has  a  repeated  frequency,  that  is, 

-  .  A  theorem  of  linear  algebra  states  that  if 

the  eigenvalue  is  repeated  "p"  times,  there  will  be  "p" 
linearly  independent  eigenvectors  associated  with  this 
repeated  eigenvalue. Since  it  is  necessary  that  these 
eigenvectors  be  orthogonal  to  each  other,  it  is  possible  to 
choose  the  eigenvectors  such  that  they  will,  in  fact, 
satisfy  the  orthogonality  relationships  of  Eqns.  (154)  and 
(155),  even  though 
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Tetble  2 

OTSTRUOL  Generated  Precpiencies  and  Nodal  Patterns 
with  Calculated  Modal  Mass  for  Composite  Graphite  Square  Plate 


Young's  Modulus=40xl0“lb/in^,  Poisson's  ratio=0.3,  Density=0 . 0542  Ib/in’ 

#-Mode  Number,  NOE-Number  of  Elements,  Final-Final  Convergent  Value,  NF-Nodal  Pattern 
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3.5  -  Coitiparison  of  the  Various  Methods  of  Determining  the 
Natural  Frequencies 

To  compare  each  method  a  tabulation  of  each  technique's 
numerical  results  for  a  free-free  thin  aluminum  plate 
appears  in  Table  1.  The  natural  frequency  values  calculated 
by  GTSTRUDL  are,  on  an  average,  3.96%  smaller  than  the 
natural  frequency  values  calculated  by  the  other  techniques. 
The  shape  patterns  obtained  by  the  methods  GTSTRUDL,  STRUDL, 
and  Warburton  were  exactly  the  same,  whereas,  Lemke ' s 
patterns  were  slightly  different.  After  this  comparison  was 
made,  it  was  decided  GTSTRUDL  would  be  used  for  the 
following  various  reasons: 

a)  computation  time  -  the  time  to  compute  the  values  is 
almost  limited  to  the  actual  time  it  takes  to  enter  the 
structure's  program  file;  it  also  controlled  by  the 
number  of  elements  in  the  system,  i.e.,  the  more 
elements,  the  larger  the  stiffness,  the  longer  it  takes 
to  compute  the  eigenvalues.  Even  with  these  limitations 
it  still  takes  GTSTRUDL  much  less  time  to  compute  the 
eigenvalues  than  with  any  other  method, 
b)  convenience  -  there  are  two  ways  in  which  the  GTSTRUDL 
method  is  more  convenient  than  the  other  methods: 
i)  The  first  one  is  that  the  actual  printed  results  are 
much  easier  to  obtain  and  use  than  the  Warburton  and 
Lemke  formulations,  i.e.,  if  it  is  necessary  to  obtain 
a  modal  deflection  at  any  particular  point,  GTSTRUDL 


readily  lists  the  normalized  eigenvector  at  each  node 
(grid  point) ;  an  interpolation  or  graphical 
representation  can  be  easily  applied  for  any  points  in 
between  the  grid  points. 

ii)  The  second  way  is  that  the  GTSTRUDL  version,  J,  is 
available  on  Howard  University's  School  of 
Engineering's  Vax  750  &  780  system;  STRUDL  is  available 
on  Howard  University's  IBM  3090  mainframe.  The  STRUDL 
that  is  available  on  the  mainframe  is  a  very  early 
version  of  the  structural  design  language;  the  version 
of  GTSTRUDL  available  on  the  VAX  is  a  more  recent  copy 
of  GTSTRUDL  than  the  STRUDL.  Also,  the  IBM  3090 
mainframe  computer  system  is  not  as  user  friendly  as 
the  Vax  750  &  780  computer  syst*  .s. 

A  listing  of  the  GTSTRUDL  commands  and  the  program 
implemented  on  the  VAX  750  appear  in  Appendix  B. 


CHAPTER  4  -  DEVELOPMENT  OF  THE  STATE  SPACE  EQUATIONS 

Feedback  control  engineering  may  be  regarded  as  the 
conscious,  intentional  use  of  the  mechanism  of  feedback  to 
control  the  behavior  of  a  dynamic  process.  This  aspect  of 
control  system  engineering  is  generally  called  control 
'theory';  whereas,  the  state-space  method  is  considered  to 
be  the  cornerstone  of  modern  control  theory. The 
advantage  of  using  state-space  methods,  instead  of  the 
frequency-domain  approach,  is  the  characterization  of  the 
processes  of  interest  by  differential  equations  instead  of 
transfer  functions. 

In  the  state-space  approach,  all  the  differential 
equations  in  the  mathematical  model  of  a  system  are  first 
order  equations;  only  the  dynamic  variables  and  their  first 
derivatives  (with  respect  to  time)  appear  in  the 
differential  equations.  The  dynamic  variables  that  appear  in 
the  system  of  first-order  equations  are  called  the  state 
variables.  The  external  inputs  are  called  the  control 
variables. 

This  chapter  will  encompass  the  development  of  the 
system's  constant  state  and  control  influence  matrices,  [A] 
and  [B],  respectively,  and  their  relationship  to  the 
system's  state  and  control  vectors,  (X)  and  {U}, 
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respectively.  Once  the  matrices  and  vectors  are  formed,  they 
will  be  put  together  to  form  the  general  state  space 
equation. 

{X}  -  [A]  m  ^  m  {U)  (156) 

In  a  linear  system  the  output  vector,  {Y),  is  assumed 
to  be  a  linear  combination  of  the  state  and  the  input, 

{y}  -  [c]  {JY}  4  [p]  (157) 

For  our  system  the  output  is  assumed  to  be  completely 
observable  and  only  a  function  of  the  state  variables,  thus, 

{y}-[C]{A'}  (158) 

where  [C]  is  an  identity  matrix. 

Before  a  complete  formulation  of  the  state-space 
equation  can  be  accomplished,  two  factors  in  the  completion 
process  must  be  discussed:  actuator  placement  and  modal 
mass. 

4.1  -  Discussion  of  the  Actuator  Placement 

When  considering  actuator  placement,  there  are  two 
points  to  consider:  (1)  orientation  control,  i.e,  providing 
the  most  effective  control  torque;  and  (2)  deformation 
control,  i.e.,  providing  the  most  efficient  actuator  force 
effort  to  control  the  shape  of  the  system. 
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Orientation  control  is  implemented  bv  the  control 
torque  vector,  Tp,*  its  mathematical  model  is  represented  by 
the  vector  cross  product  of  the  actuator  position,  vector  and 
the  actuator  force  vector.  To  acquire  the  maximum  torque  in 
any  direction  one  may  either  increase  the  force  of  each 
actuator  and/or  increase  the  torque  arm.  Increasing  the 
force  that  contributes  to  the  control  effort  may  result  in 
an  undesirable  increase  in  the  energy  required.  An  increased 
torque  arm  can  be  realized  by  placing  the  actuators  at  some 
maximum  distance  from  the  origin;  this  is  plausible  as  long 
as  these  placement  positions  do  not  fall  on  the  modal  nodal 
lines  (the  zero  deflection  lines  of  the  modal  patterns) . 

This  leads  to  the  second  consideration,  deformation 
control.  Placement  of  an  actuator  on  a  modal  nodal  rine  must 
be  avoided.  An  actuator  placed  on  a  nodal  line  of  a 
particular  mode  will  have  no  effect  in  contributing  to  the 
defoirmation  control  of  that  mode.  Noting  that  each  mode  has 
a  nodal  line  pattern,  it  is  important  to  accurately 
determine  each  mode's  pattern  shape  (see  Table  2). 

The  nodal  lines  of  the  first  two  modes  for  a  thin 
square  plate  are  easy  to  locate;  unfortunately,  the  third 
mode's  nodal  circle  is  more  difficult.  With  the  use  of 
GTSTRUDL-144  BPR  element  data,  the  third  mode's  eigenvectors 
-for  the  nodes  located  on  the  midway  line  of  the  plate's 


major  surface-  vvere  plotted.  By  viewing  the  plot,  vsee 
Figure  5),  as  a  cross-sectional  cut  midway  through  the 
plate,  the  determination  of  the  nodal  circle  associated  with 
the  third  mode  can  be  accurately  located. 


Third  Mode  Normalized  Deflection 
For  Square  Plate  At  Midline 


Normalized  Deflection 

0.2  I - — - - 
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Figure  5.  Normalised  Deflection  Plot  for  3^^  Mode's  Cross- 
Sectional  View 


Once  the  nodal  line  pattern's  were  acquired  for  the 
first  three  inode's,  a  second  plot,  (see  Figure  6),  was  drav.'n 
indicating  3  possible  placements  of  the  actuators  (i.e., 
etc.)*  Note,  in  previous  references  [18,19]  two 
different  sets  of  actuator  positions  were  studied:  these 
proved  to  be  less  effective  because  of  their  placement  on 
nodal  lines  (see  Figure  7).  Two  sets  of  actuator  positions, 

A  &  A'  have  been  chosen  for  this  thesis  (see  Figure  8).  Six 
ac-uators  are  assumed  to  be  placed  in  a  symmetrical  pattern 
with  both  force  axis  directions  perpendicular  to  the  major 
surface  and  along  the  edge  of  the  plate.  The  first  four 
actuators  are  assumed  to  be  placed  normal  to  the  major 
surface;  they  help  control  the  shape  deformation  and  the 
torque  about  the  two  major  axes  in  the  plane  of  the  plate. 
The  other  two  actuators,  5  and  6,  are  placed  along  the  edge 
of  the  plate;  they  provide  control  for  torque  about  the 
third  axis,  normal  to  the  major  surface. 


Figure  6.  144  'BPR'  Finite  Element  Grid  Pattern  with  Nodal 

De-flections  for  Three  Different  Sets  of  Actuator  Positions 


Figure  7.  Previously  Studied  Actuator  Locations,  sets  I  4  II 


To  calculate  the  modal  mass  values,  a  formulation 
presented  by  Hurty  and  Rubenstien  was  followed  - Using  the 
equation  of  the  form: 

n. o. e. 

n-l 

where  r  represents  the  mode  number;  n  represents  the 
gridpoint  node;  n.o.e.  represents  the  total  number  of 
elements;  M,  represents  the  mass;  and  ^  represents  the 
normalized  deflection  of  the  grid  point  (node) . 

The  modal  masses  were  calculated  and  plotted  for  the 
plate  modeled  by  different  numbers  of  assumed  finite 
elements.  Figure  9  shows  the  variation  of  the  modal  mass  as 
a  function  of  the  number  of  finite  elements  assumed  to 
represent  the  system.  The  final  convergent  value  was  used 
for  the  nodal  mass.  The  calculated  modal  masses  for  the 
first  three  modes  are,  M, =20, 278. 65  kg,  M2=29,366,14  kg,  and 
M3«=18572.3  4  kg. 


Calculated  Modal  Mass  -1st,  2nd,  &  3rd- 
versus  Number  of  Elements 


Modal  Mass,  Mr  (slugs) 


1st  Modal  Mass  +  2nd  Modal  Maas  3rd  Modal  Maas 


Figure  9.  Modal  Mass  versus  Number  of  Finite  Elements 
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4.3  -  Formulation  of  the  System's  State  Vector  and  Matrix 


With  the  assumption  that  the  system  can  be  modeled  by 
three  rigid-body  rotational  modes  and  the  first  three 
generic  modes,  assignment  of  the  correct  dynamic  variables 
to  the  state  vector  follows.  The  state  vector  for  this 
system : 


where  the  state  variables  are  defined  as, 

4>  -  X,  4)'  -  -  X,  -  x, 

^  at  ^ 


ij;  -Xj 

V-Xg 

e  -  Xj 

B'-Xg 

Zi  -x^ 

Zi-^xo 

Z2-X5 

Z2  “ 

^3-^6 

Zx  “  X.  2 

(160) 


(161) 


Preparation  of  the  system's  state  matrix  through  the  usage 
of  Eqns.  (82-85)  follows. 


lA]  - 


-4Q 

0 

0 


0 

0 


to], 

0 

0 

3C, 


3  - 


0 

0 

0 


-fJKiv 


0  0  0 
0  0 
0  0 

0  0 


(162) 
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With  the  substitution  of  the  parameters  for  the  nominal 
orientation  of  Case  (1),  (see  Figure  2.1);  the  principal 
moment  of  inertias  I, -2. 5x10®  kgm^,  1=1=1.25x10®  kgm^;  the 
angular  velocity  <i)j.=0 . 0011162  rad/s;  and  the  modal 
frequencies  a, =0.0547  rad/s,  W2=0. 07852  rad/s,  and  03=0.09773 
rad/s,  the  system's  state  matrix  [A],  becomes: 

"  •  f  6x6  •  •  •  ■  [-T]  6x6 

400  0  0  00  -2  0 

000  0  0  010... 

[A]  -  0  0  3  0  0  0  0 

00  0  -2398  0  0:0. 

0  0  0  0  -4945  0  :  . 

000  0  0  -7663  0  0 

4.4  -  Formulation  of  the  System's  Control  Vector  and  Matrix 

The  system's  control  effort,  [B]{U),  is  defined  by  the 
matrix  multiplication  of  the  system  control  matrix,  [B],  and 
the  system  control  (input)  vector,  (U).  This  term  in  the 
state-space  equation  allows  for  control  feedback  for  the 
system.  The  general  control  influence  matrix,  [B] ,  is 
defined  as, 

^  ®  ^  6x6 

’l[^]6x6. 


(164) 


where  the  lower  part  of  the  [B]  matrix  depends  on  the 
actuator  positions.  The  control  effort,  [B]{U),  is 
formulated  as  such. 


[B]  {U]  - 


[  Tp/ti)  cl  -ix6 


(165) 


where  Tp  represents  the  external  acceleration  due  to  the 
actuator;  and  represents  the  external  force  on  the  n— 
mode,  their  formulation  is  as  follows: 


(16b) 


-^p  ^p 

1-^*3 


(167) 


where  Npj  represents  the  control  torque  on  the  location 
due  to  the  j—  actuator;  fp  represents  the  position  vector 
from  the  origin  to  the  p^  actuator;  and  f^  represents  the 
control  force  vector  due  to  the  j—  actuator.  For  the 
generic  modal  equations,  the  control  forces  can  be 
transformed  into  the  corresponding  modal  forces 


E,  -  f ^  (z^) an?  (168) 

J  P 

where  represents  the  modal  shape  function  at  the 

location  corresponding  to  the  n—  mode. 


To  define  the  system  control  matrix  for  the  nominal 
orientation  of  Case  (1),  it  is  necessary  to  establish  the 
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correct  position  and  force  vectors  for  this  system.  For  the 
actuator  sets,  A  and  A',  the  general  position  vectors  and 
force  vectors  are: 

jf,  .  -  yj  -  zJc  :  f ,  .  - 

1-4  1-4 

-^5, 6 

The  angular  acceleration  produced  by  each  actuator  becomes: 

/  ly- yi/ 1  j,)  (170) 

The  lower  half  of  the  system  control  matrix  contains 
external  force  terms.  A  definition  of  the  system's  modal 
deflections  at  each  actuator  location  is  necessary,  along 
with  the  modal  mass.  Because  the  actuators  are  placed  in  a 
symmetrical  pattern,  the  deflections  for  actuators  1-4  are 
equal  in  magnitude  but  alternate  in  sign.  Since  actuators  5 
and  6  are  located  on  the  side  of  the  plate  their  deflections 
are  assumed  to  be  zero.  Thus  the  control  force  terms  become; 


The  system's  control  effort,  for  both  sets  of 

actuator  locations,  A  and  A',  follows: 
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[B]  {f/} 


0 

0 

0 

0 

-  .  1605 

.1605' 

'A 

£ 

- . 1605 

.321 

.  1605 

-.321 

0 

0 

^2 

-  .  321 

-  .  1605 

.  321 

.  1605 

0 

0 

-  .2221 

.  2221 

-.2221 

.2221 

0 

0 

-  .  189 

.  189 

-  .  189 

.189 

0 

0 

■fs 

.1842 

.  1842 

.  1842 

.1842 

0 

0 

A 

[B]  {U]^> 


0  0 
.1605  .214 

-.214  -.1605 

.1545  .1545 

.0555  .0555 

0306  .0306 


0  0 
.1605  -.214 

.214  .1605 

-.1545  .1545 
-.0555  .0555 
.0306  .0306 


-.1605  .1605' 
0  0 
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0  0 

0  0 

0  0 


A 

A 

f 


i  (172) 


)  (173) 


CHAPTER  5  -  PRESENTATION  OF  LINEAR  QUADRATIC  OPTIMUM  CONTROL 
THEORY 


The  main  purpose  of  linear  control  theory  is  to  design 
the  correct  compensator  (gains)  for  the  given  system.  For 
this  thesis  an  application  of  optimum  control  theory  will  be 
implemented  to  synthesize  the  control  law  gains. 

In  this  chapter,  as  an  optimum  control  theory  algorithm 
the  linear  quadratic  regulator  (LQR)  technique  was  chosen. 
Initially,  the  performance  criteria  is  discussed,  the 
system's  linear  differential  equation  is  solved,  and  the  LQR 
technique  is  stated  for  a  continuous-time  system.  Finally, 
the  sampling  period  criteria  for  the  discretization  of  a 
continuous-time  control  system  is  presented;  following  the 
sampling  period  criteria,  an  application  of  the  LQR 
technique  to  a  discrete-time  system  is  performed. 

5.1  -  Discussion  of  the  Performance  Criteria  or  Cost 
Function 

The  dynamic  process  considered  here  is  characterized  by 
the  vector-matrix  differential  equation 

x(  t)  -  A  ( t)  x  (  t)  4  B  (  t)  U  (  t)  (174) 
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where  the  variables  have  been  previously  defined  in  Chapter 
4.  This  chapter  will  seek  to  define  the  linear  control  law 

u{t) x(t)  (175) 

where  is  a  suitable  gain  matrix  related  to  u(t).  An 
attempt  will  be  made  to  find  the  optimum  gain  matrix  to 
minimize  the  performance  criterion,  J,  (or  "cost  function") 
expressed  as  the  integral  of  the  quadratic  form  in  the  state 
error,  e(t),  plus  a  second  quadratic  form  in  the  control, 
u(t) ;  i.e. , 

,  Fe(tf)  >  +  ^  r<e(t)  ,  Q{  C)  e  { t)  >  +  <i:  { t)  ,  Rit)  ui  C)  >] 


Initially,  there  are  some  assumptions  to  be  made.^^’^ 

a)  The  terminal  time,  T,,  is  specified. 

b)  F  is  a  constant  mxm  positive  semidef inite  matrix. 

c)  Q(t)  is  an  mxm  positive  semidef inite  matrix. 

d)  R(t)  is  an  rxr  positive  definite  matrix. 

Previously  mentioned  was  the  assumption  that  the  system  is 
completely  observable,  thus  C{t)=I.  The  error  vector  is 
defined  by 


e  ( t)  -  z  ( t)  -  y  ( t) 


(177) 
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where  z(t)  is  assumed  to  be  zero.  Therefore, 

-e  i  t) -yi  t)  -  x(  C)  (178) 

Now  the  control  function,  J,,  can  be  designated  as, 

J^~^<xiTf)  [<x{t)  ,  0(c)  x{t)>^<ui  c)  ,R(C)u(t)>]dt(X79) 

Before  an  attempt  to  find  the  optimum  gain  matrix,  is 

made,  some  comments  about  the  oost  function  are  in  order. 
With  respect  to  the  limits  on  the  integral,  the  lower  limit, 
tjj,  is  identified  as  the  present  time,  and  the  upper  limit, 
Tf,  is  the  terminal  time  or  final  time.  The  time  difference, 
T^-tg,  is  the  control  interval,  or  ' time-to-go ' ,  If  the 
terminal,  T^  is  finite  and  fixed,  the  time-to-go  heeps 
decreasing  to  zero,  at  which  time  the  control  process  ends. 
However,  in  the  customary  case,  the  terminal  time  is 
infinite.  In  this  case  we  are  interested  in  the  behavior  of 
the  process  'from  now  on',  including  the  state. 

Attention  is  now  focused  on  the  term; 
l/2<x  (T^)  ,  Fx  (T^)  > .  This  term  is  often  called  the  terminal 
cost;  its  purpose  is  to  guarantee  that  the  final  state  is 
small  at  the  terminal  time.  This  should  be  included  in  the 
final  state,  if  x(T^)  is  expected  to  be  particularly  large. 
Otherwise,  F  can  be  set  to  zero  and  the  rest  of  the  cost 
function  can  be  relied  upon  to  guarantee  that  the  terminal 
state  is  not  excessively  large. 
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Finally,  consider  the  weighting  matrices,  Q  and  R. 

These  are  often  called  the  state  weighting  matrix  and  the 
control  weighting  matrix,  respectively.  A  formula  fc»-  the 
control  gain  matrix,  G^,  can  be  derived  such  that  it 
involves  terms  of  the  weighting  matrices  and  the  steady 
solution  to  the  Ricatti  matrix  eqpaation.  By  plugging  in  the 
matrices  Q  and  R  -along  with  the  matrices,  A  and  B,  that 
define  the  dynamic  process-  into  the  computer  routine, 
ORACLS,  can  be  found.  If  the  process  is  controllable 
and  Q  and  R  are  suitable,  the  computer  finds  a  suitable  G^. 
This  is  not  to  say  that  the  calculation  is  a  trivial  problem 
-far  from  it-  but  only  that  the  problem  of  determining  G^ 
once  A,  B,  Q,  and  R  are  given,  is  not  a  control  design 
problem  but  a  problem  in  numerical  analysis. 

The  question  of  concern  to  the  control  system  designer 
is  the  selection  of  the  weighting  matrices,  Q  and  R.  In  the 
cost  function  defined  by  Eqn.  (179),  two  terms  contribute  to 
the  integrated  cost  of  the  control:  the  quadratic  form  x^Qx, 
represents  a  penalty  on  the  deviation  of  the  state  x,  from 
the  origin  (this  means  that  the  desired  state  is  at  the 
origin,  F=0,  not  at  some  other  state)  and  the  term  u^Ru, 
represents  the  'cost  of  control'.  The  physical 
interpretation  of  is  this:  We  wish  to  keep  the  state  near 
zero  without  excessive  control -energy  expenditi’.re. 


85 


The  weighting  matrix,  Q,  specifies  the  importance  of 
various  components  of  the  state  vector  relative  to  each 
other.  For  example,  suppose  x,  represents  the  system  error, 
and  represents  the  successive  derivatives,  i.e.,'^^^ 

x^-x 

Xj  -  X 

X3-X  (180) 

Xy~X>^-^ 


If  the  error  and  none  of  its  derivatives  are  of  concern, 
then  one  might  select  a  state  weighting  matrix  such  as: 

[1  0  •  0 
0  0-0 

0  0-0 

which  will  yield  in  the  quadratic  form 

x'^qx-xl  (182) 

But  the  choice  of  Eqn.  (181)  as  a  state  weighting  matrix  may 
lead  to  a  control  system  in  which  the  velocity  X2=x  is 
larger  than  desired.  To  limit  the  velocity,  the  performance 
integral  might  include  a  velocity  penalty,  i.e., 

X  ^gx- xl  +  c^xl  (183) 

which  would  result  from  a  state  weighting  matrix 


(181) 


0 


1  0  ■■  ol 

I 

0  C“  •  Oi 
0  0  ■  0; 


(184) 
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The  choice  of  the  state  weighting  matrix,  Q,  depends  on  what 
the  system  designer  is  trying  to  achieve.  The  considerations 
of  the  above  with  regard  to  Q  apply  to  the  control  weighting 
matrix,  R.  The  term,  u’Ru,  in  the  perfonnance  index,  J,,  is 
included  in  an  attempt  to  limit  the  magnitude  of  the  control 
signal,  u(t).  For  this  thesis,  Q  and  R  are  chosen  to  be 
constant  matrices,  where 

The  relationships  between  the  weighting  matrices,  Q  and 
R,  and  the  dynamic  behavior  of  the  closed  loop  system  depend 
on  the  matrices  A  and  B  and  are  quite  complex.  It  is 
difficult  to  predict  the  effect  of  a  given  pair  of  weighting 
matrices  on  the  system’s  closed-loop  behavior.  A  suitable 
approach  for  the  designer  would  be  to  solve  for  the  gain 
matrices  that  result  from  a  range  of  the  weighting  matrices, 
Q  and  R,  and  calculate  (or  simulate)  the  corresponding 
closed  loop  response.  The  gain  matrix,  G^,  that  produces  the 
response  closest  to  meeting  the  design  objectives  is  the 
ultimate  selection.  With  the  ORACLS  software,  it  is  a  simple 
matter  to  solve  for  given  A,  B,  Q,  and  R.  In  a  few  hours 
time,  the  gain  matrices  and  transient  response  that  result 
for  a  dozen  or  more  combinations  of  Q  and  R  can  be 
determined,  and  a  suitable  selection  of  G^^  can  be  made. 
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5.1.1  -  Solution  of  the  Linear  Differential  Equations  in 
State-Space  Form 


Starting  with  the  system  differential  state  equation, 


xit)  -Axi  t)  Bu(t) 


where  the  matrices  A  and  B,  are  considered  constant 


(187) 


matrices,  the  simplest  form  of  the  general  differential 
equation,  Eqn.  (186)  is  the  homogeneous,  i.e.,  unforced 


equation 


X { t)  ~Ax{  t) 


The  solution  can  be  expressed  as 


x(t)  - ce 


(188) 


(189) 


where  c  is  a  suitably  chosen  constant  vector.  To  verify  Eqn. 
(189)  the  derivative  is  calculated  as 


•cAe'^'  -Ax{t) 
at 


(190) 


To  evaluate  the  constant,  c,  suppose  that  at  some  time,  r, 
the  state  x(t),  is  given  as 

x{x)  -  (191) 

Multiplying  both  sides  by  the  inverse  of  e*’^  leads  to 


c-  (e'^'')  "^x(t)  (192 ) 


Substitution  into  the  homogeneous  solution  leads  to 
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x(  t)  -  (193) 

where  the  matrix  is  a  special  form  of  the  state- 

transition  matrix. 

Through  the  use  of  'method  of  variation  of  the 
constant',  a  particular  solution  to  the  nonhomogeneous ,  or 
'forced'  differentia]  equation  is  found.  Seeking  a  solution 
of  the  form  of  Eqn.  (187)  one  can  select 

x(  C)  -  e  t)  (194  ) 

where  c(t)  is  a  function  of  time  to  be  determined.  After 
taking  the  time  derivative  of  x(t)  and  substituting  into 
Eqn.  (185)  the  following  is  obtained 

Ae^’^c(t)  *  (195) 

Upon  cancellation  of  the  terms  Ae*'c(t)  and  multiplication 
of  the  remainder  by  e'**, 

d{t)  -  (196) 

Thus,  the  desired  function  c(t),  can  be  obtained  by  simple 
integration 

£ 

c  (  c )  -  je  '^Bu  ( A. )  dl 

T 


(197) 
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The  lower  limit,  T,  on  this  integral  is  undefined  for  now; 
substitution  into  the  particular  solution  produces 

c  t 

x(t)  -  e  e'^^Buik)  dX-  Bu  f  k)  dk  (198) 

r  T 

The  complete  solution  to  Eqn.  (187)  is  obtained  by 
adding  the  'complementary  solution',  Eqn  (193),  to  the 
particular  solution,  Eqn.  (198).  The  result  is 

t 

x(  c)  -  x{x)  Bu  ik)  dk  (199) 

T 

The  proper  value  for  the  lower  limit,  T,  on  the  integral  can 
now  be  determined.  At  t=T,  the  complete  solution  becomes 

t 

x( t) -x(t)  *  Je'''^‘*’Bu(A.)  dl  (200) 

T 

The  integral  in  Eqn.  (200)  must  be  zero  for  any  u(t);  this 
is  only  possible  if  T=t .  Therefore,  the  complete  solution  to 
Eqn.  (187) ,  when  A  and  B  are  constant  matrices, 

C 

x{t)  -e^‘‘'*’x(T)  *  Bu(k)  dk  (201) 


5.2  -  Discretization  of  a  Continuous-Time  Control  System 


When  a  continuous-time  control  system  with  complex 
poles  is  discretized,  the  introduction  of  sampling  may 
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impair  the  controllability  and  observability  of  the 
resulting  discretized  system.  That  is,  pole-zero 
cancellation  may  take  place  in  passing  from  the  continuous¬ 
time  case  to  the  discrete-time  case.  Thus,  the  discretized 
system  may  lose  controllability  and  observability. 

A  system  which  is  completely  state  controllable  and 
observable  in  the  absence  of  sampling  remains  completely 
state  controllable  and  observable  after  the  introduction  of 
sampling,  if  and  only  if,  for  every  eigenvalue  of  the 
characteristic  equation,  the  relation 

Re  -  Re  Xj  (202) 

implies 

Im  iX^-Xj)  *  (203) 

where  are  the  eigenvalues  oj.  t.he  continuous-time  system 
matrix  A,*  T^  is  the  sampling  period,  and  n=±l ,  ±2 ,  .  .  . 

In  addition  to  the  f orementioned  requirement,  two 
points  should  be  taken  into  consideration  when  choosing  a 
sampling  period 

a)  If  the  sampling  interval  is  too  long,  the  performance 
of  the  sampled  data  system,  tends  to  deteriorate;  this 
makes  signal  reconstruction  difficult. 

b)  Implementation  of  a  very  short  sampling  interval  may  be 
limited  by  computer  operation  times  and  the  expense  of 
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fast  A/D  and  D/A  converter  devices;  thus,  overtaxation 
of  the  computer  system  with  data  may  result. 

Therefore,  the  sampling  time  should  be  as  large  as  possible 
after  the  performance  of  a  sampled  data  system  meets  the 
requirements  of  the  design.  Using  the  open-loop  eigenvalues 
calculated  by  ORACLS  a  table  was  developed  listing 
unacceptable  sampling  periods  (see  Table  3  and  4) . 


Table  3 

Eigenvalues  and  Moduli  for  Continuous-Time  Open-Loop  System 
Calculated  by  ORACLS  for  Orientation  of  Case  (l) 


Continuous  Open  Loop 

Eigenvalues 

Moduli 

Real 

Imaginary 

■1 

0 

0 

0 

2 

0 

0 

0 

3 

1.4142 

0 

1.4142 

■1 

-1.4142 

0 

1.4142 

5 

1.7321 

0 

1.7321 

6 

-1.7321 

0 

1.7321 

IB 

0 

49.006 

49.006 

8 

0 

-49.006 

49.006 

9 

0 

70.346 

70.346 

10 

0 

-70.346 

70.346 

11 

0 

87.556 

87.556 

12 

0 

-87.556 

87 . 556 

Ci  n 


Table  4 

Calculation  of  unacceptable  Szunpling  Periods 
from  the  System  Matrix  Eigenvalues 


i\j 

0 

-49.00 

♦49.00 

-70.35 

♦70.35 

-87.55 

♦87. SS 

0 

0 

114.8 

-114.8 

-80.00 

64.23 

-64.23 

-49.00 

-114.8 

0 

-57.42 

263.7 

-47.12 

145.9 

-41.21 

+49.00 

114.8 

57.42 

0 

47.12 

-263.7 

41 .21 

-145.9 

-70.35 

-80.00 

-263.7 

-47.12 

0 

-39.95 

327.4 

-35.65 

♦70.35 

80.00 

47.12 

263.7 

•39.95 

0 

35.65 

-327.4 

-87.55 

-64.23 

-145.9 

-41.21 

-327.4 

-35.65 

0 

-32.07 

♦87.55 

64.23 

41.21 

145.9 

35.65 

327.4 

32.07 

C 

all  nuTbers  a'-e  multiplied  by  It, 

k=^+/-1,2,3.. 

i  -  Ith  Imaginary  Eigenvalue 
j  -  Jth  Imaginary  Eigenvalue 


5.2.1  -  The  Linear  Quadratic  Regulator  Technicrue  Applied  to 
a  Continuous-Time  System 

Let’s  begin  by  considering  the  linear  time-varying 
system  of  Eqn.  (202), 

x{t)  -Ait)  xit)  Bit)  uit)  (204) 

the  cost  function,  J^,  given  by  Eqn.  (205), 

<xi  t  f)  ,  Fxi  c, J"  I  <x(  t)  ,  t)  x(  t'l  ( t)  ,  B  (  C)  u{  t)  >]  dt  (205) 


and  the  optimal  control  u{t),  is  selected  such  that  the  cost 


function  is  minimized . It  is  defined  as, 


(206) 


u(t)  -  -i?'Mc)  K(t)  x(  C) 

where  K(t)  is  a  nxn  symmetric  matrix  and  is  the  unique 
solution  of  the  Ricatti  equation^’’ 


--KiOAiC)  -A  ^  it)  Kit)  *  Kit)  Bit)  R-^  it)  B^it)  Kit)  -Qit)  (207) 


which  satisfies  the  boundary  condition 


K(Tf)  -F 


(208) 


The  state  of  the  optimal  system  is  then  the  solution  of  the 


linear  differential  equation 


x(t)  -  [A(t)  -  ]  x(c) 


(209) 


The  block-diagram  representing  the  optimal  control  system  is 
displayed  in  Figure  10. 


Figure  10.  Block  Transfer  Function  of  the  Optimal  Control 
System 
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The  response,  x(t)  ,  of  the  optimal  system  can  now  be  w^ritten 
as  the  solution  to  the  differential  equation 

x{t) -Git)  x(t)  (210) 

where  G(t)  is  an  nxn  matrix  given  by 

Git)  -  Ait)  -  Bit)  R-^  it)  it)  Kit)  (211) 

To  find  the  gain  matrix,  K(t) ,  a  steady  state  solution 
of  the  matrix  Ricatti  equation  (Eqn.  (207))  is  necessary, 
when  the  final  time,  T^,  approaches  infinity;  thus,  the 
derivative  K(t)=o,  making  K (t) ^Constant  (see  boundary 
conditions,  Eqn.  (206)).  Since  the  Ricatti  equation  is 
nonlinear,  closed-form  solutions  usually  cannot  be  obtained; 
therefore,  K(t)  must  be  computed  using  a  digital  computer. 
For  the  purposes  of  this  thesis  the  ORACLE  subroutines  will 
be  implemented. 

ORACLE  is  a  modern  ccntrol  theory  design  package  for 
constructing  controllers  and  optimal  filters  for  systems 
modeled  by  linear  time-invariant  differential  or  difference 
equations.  The  digital  FORTRAN-coded  ORACLE  system 
represents  an  application  of  some  of  today's  best  numerical 
linear  algebra  procedures  to  implement  the  Linear-Quadratic- 
Regulator  or  Guassian  methodology  of  modern  control  theory. 
An  example  of  the  simulation  program  which  incorporates 
various  ORACLE  subroutines  has  been  placed  in  Appendix  D. 
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One  important  thing  to  realize  is  that  the  gain  matrix^ 
K(t),  is  independent  of  the  state,  so  that  once  the  system 
and  the  cost,  J,  ,  have  been  specified,  K(t)  can  be  computed 
before  the  optimal  system  starts  to  operate. 

5.2.2  -  Formulation  of  the  Linear  Quadratic  Regulator 
Technique  Subject  to  a  Discrete-Time  Domain 


Consider  the  continuous  time  control  system 


X  it)  -  Ax  { t)  Bu  (  c) 


subject  to 


(212) 


r 

-ix  ’’{ tf)  Fx(  Cf)*~f[x'‘'(c)C’x(C}-i-u^(c)Bu(t:)]dt  (213) 

0 

If  the  continuous-time  control  system  is  approximated  by  its 
discrete  equivalent,  the  sampling  period  T^,  can  be 
represented  as 

~~kT^,  k-0,1.2,  .,N  (214) 

thus,  its  discrete  equivalent  becomes. 


x(  ik^l)  T^)  -GiT^)  x(kT^)  uikT^)  (215) 

and  the  discretized  performance  index,  when  the  final  time 
T,=NT^,  becomes 


K-l 


J-  Ix^iNT^)  FxiNT^)  O^xikr,) 

^  ^  *r-0 

+  ?.x  '^(kT^)  M^u(kT^)  ■*-  u  '^(kTj  R.  u(kT^)  ] 


(216) 


It  is  noted  that  the  integral  term  in  Eqn.  (213)  is  net 
repJ  aced  by 


—  [x^ikl'j  QxikT^)  -u 

“  *-o 


(217) 


but  is  modified  to  incluae  a  cross  term,  involving  x(kT^)  and 
u(k:T^).  Also,  matrices  Q  and  R  are  modified.  By  considering 
the  discretized  quadratic  optimal  control  problem  by  use  of 
a  simple  example  -similar  to  the  systeTt;  considered  in 
section  5.1.1-  the  terms  Q, ,  K, ,  and  R,  are  developed. 


Consider  the  continuous-time  system  defined  by 

X  i  t)  •  Ax  i  t)  ■*  Bu  K  L)  (218) 

where  A  and  B  are  constant  matrices  and 

'J(t)  -uikT^)  .  ik*l)  (219) 

The  performance  index  to  be  minimized  is 

J-^ -jxUnt^)  ^  ^  J  [  CxM  t)  -  R’u’M  t)  ]  dt  (220) 


To  begin,  the  sv'stero  equation  and  performance  index 
must  be  transformed  to  a  discrete-time  domain;  afterwards 
the  discretized  quadratic  optimal  control  problem  can  be 
formulated.  Equation  (218)  may  be  discretized  as  follows: 


x(  {k*l)  TJ  -GiT^^)  xikT^)  +/V(rj  uikT^) 


(221) 


where , 


G{T,) 


Ar, 


7*. 

HiT^)  -  j 


or 


x{  ik^l)  TJ  ~e^^‘x(kTJ  +  —  uikT) 

s  s  ^  s 

The  performance  index,  J,  given  by  Eqn.  (220)  may  be 
discretized.  First  rewrite  J  as 


s-i 

J- ^xUNT^}  ^  f  (C’XM't)  +FuMt)]dt 


jtr. 


Noting  that  the  solution  x{t)  for  kTg<t<  (k+l) can  be 
written  as 


x(t)  ♦  I  e^‘^-’>Su{T) 

*r. 


dr 


-5  { t-kTJx{kTJ  +r|  ( t-kr.)  L'(kr-) 


where , 


iit-kTJ  -e 


A(  t->cT.) 


r|  ( t-kTj  “  /  ^  ^  " 

kT, 


B\^AU-kT, 


-  l' 


(222) 


(223) 


(224) 


(225) 


(226) 


(227) 


(226) 


The  perfommance  index  J,  can  be  written  as  follows 


W-l  ‘ 

*-c  jj 

+  11  (  C-kTj  uikT^)]^*Ru'UkT^)]dt 

-^x^(NT^)  -~Y,  l(^,x-ikT^)  *2M,x{kT^)  uikT^)  *R,u'-(kTJ] 

2  ^  k-o  ' 


(229) 


where, 


'.k.  11  r. 


Pj-  J  C>$2(t-/cr^)dt 


ik-V.Z, 


fj-  J  Qi(t-kT^)r)  {t-kT,)  cit 


(k-l)T, 

i?i-  f  [C>r\^(t-kT,)  *R]dc 


(230) 


(231) 


(232) 


Notice  that  Q,,  M, ,  and  R,  may  be  simplified  as  follows: 

c.-  /  ( 

kT. 


(233) 


(k*l)  T, 


1,-  /  Pe"”"- 


A!  r  kT,  R  r  _y5l 


•  ’  -lidt-  1) " 

2A^ 


(234) 


<*‘i)  r, 

fl.-  I  [p{^u'"-"-'-i])‘ 


‘  J  riA 

kT,  ‘ 

D2r)r  .»• 


(235) 


1)  +  2,A7J  +  /?T. 
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Summarizing,  the  present  discretized  quadratic  optimal 
control  problem  may  be  stated  as  follows.  Given  the 
discretized  system  equation 

x{  (k^l)  T^)  ‘G{T^)  x(kT^)  u(kT^)  (236) 

where , 


and  H(T^)  - -f  (237) 

find  the  optimal  control  sequence,  u(0),  u  (T^)  , . . . , u ( (N- 
1)T5)  such  that  the  following  performance  index  is  minimized 


N-l 


j-  -1x2  (JcTj)  *2M^x{kT^)  u(kTg)  ikT^)  1 


(238) 


k-O 


Such  a  performance  index  including  a  cross  term  involving 
x(kTg)  and  u(lcTg)  can  be  modified  to  a  form  that  does  not 
include  a  cross  term,  and  the  solution  to  the  discretized 
quadratic  optimal  control  problem  can  then  be  obtained.  This 
subject  is  presented  in  the  following. 


Ta)cing  into  consideration  the  quadratic  optimal  control 
problem,  where  the  system  is  given  by 

x(k-^l) -Gxik)  +Hu{k)  ,  x(0)-c  (239) 

and  the  performance  index  is  given  by 

J-  —x^iN)  FxiN)  +  —  Y'  [x^(k)  O^xik)  +2x^(k)  M^u(k)  *  u  (k)  R^u  (k)  ] 
2  2 


(240) 


IOC- 

where  Q,  and  F  are  nxn  positive  definite  or  semidefinite 
Hermitian  matrices,  R,  is  an  rxr  positive  definite  Hermitian 
matrix,  and  M  is  an  nxr  matrix  such  that  the  matrix 


0^  M, 

is  positive  definite.  This  means  that 


(241) 


. .  ' 

-X  ^(ic)  O^x(k)  X  ^(k)  M^u  (k)  +  u  ’’(Jc)  M-Jxik)  *  u  ^(k)  R^u  ik) 
~x^(k)  C>^x(k)  *2x^(k)M^u{k}  *  u'^  {k)  R.^u{k) 


(242) 


is  positive  definite.  Note  that  the  performance  index,  J,  of 
Eqn.  (240)  includes  a  cross-term  similar  to  Eqn.  (242) . 


In  order  to  obtain  the  optimal  control  vector  u(k),  let 
us  define 

(243) 

and  eliminate  Q  from  the  performance  index  J,  Then  Eqn. (240) 
becomes 


J-  ^x’^m  Fxm  {X^ik)  x(jc) 

+  2x’'(ic)Wj,  u(Jc)  +  u"^(k)  R,  u(k) ) 

-^x'^(N)Fxm  [x’'(/c)£?x(Jc)  *x'^{k)M^R{^M^Jxik) 

*2x'^ ik)  M^u{k)  *  u  ^(k)  R,^u(k)  ] 

^-1 

-  ■^x’^(N)  Fx(N)  *~Y. 

^  ^  Jc-0 

+  [x  "^{k)  M,^R{^  *  u  ~ik}  ]  [R^^MJx(k)  *  u(k)  ]  } 
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Define 

v{k)  •  Ri'^’ M-Jxik)  ■*  u(k)  (245) 

such  that 

N-l 

J- iN)  Fx  (N)  —'y'  [x'^ik)  Ox(.k)  v'^ik)  R,v{k)  ]  (246) 

^  *-0 


Notice  that  Eqn.  (246)  no  longer  involves  the  cross  term;  it 
has  been  effectively  eliminated. 


By  substituting  Eqn.  (245)  into  the  system  equation, 

Eqn .  (239)  becomes 

x(A:+l)  -Gxik)  *H[v(k)  -  ] 

•  iG-HRi^M/)  Xik)  *Hv(k)  {2A1) 

-6x{k)  ■^Hv(k) 

where , 

G-G-HRl'^M^  (248) 

Note  that  the  quadratic  optimal  control  of  the  system  given 
by  Eqn.  (239)  with  the  performance  index  given  by  Eqn.  (240) 
is  equivalent  to  the  quadratic  optimal  control  of  the  system 
given  by  Eqn.  (247)  with  the  performance  index  given  by 
Eqn. (246).  Hence,  the  optimal  control  vector  v{k) ,  that 
minimizes  the  performance  index  given  by  Eqn.  (246)  can  be 
given  as  follows.  Defining 

v(k)  --  -  [R^^  P(k*l)  H]-'^  H'^P(k*l)  Gx{k)  (249) 


10 


where  P(k)  is  a  modified  version  of  the  Ricatti  Equation. 


Pik)  -O*  G'^F(k*l)  "P(k*l)  r'G  (250) 

PIN)  -F  (251) 

The  optimal  control  vector,  u(k) ,  can  then  be  given  by 

u(k)  -vik)  -Fr^Jxlk)  (252) 

where  v(k)  is  given  in  Eqn.  (249) ,  Eqn. (252)  can  be  reduced 
to  the  following  form: 

u(k)  H]  [H^Pik*!)  G->-M./]x(k}  (253) 


5.3  -  Solution  bv  the  Conventional  Minimization  Method  Using 
Lagrange  Multiplier 

In  the  present  optimization  problem,  the  minimization 
of  J  given  by  Eqn.  (246),  repeated  here 

J--  -~X  (N)  Fx{N)  ’'(^^  0x(k)  +  V  ^(k)  Rj  v(Jc)  ]  (2  54  ) 

when  it  is  subject  to  the  constraint  equation  specified  by 
Eqn.  (247) 


x(k-^l)  -(Sxik)  *Hv(k) 


(255) 
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with  the  initial  condition  on  the  state  vector 


>r{0)  -  c  (256) 

is  considered  by  using  a  set  of  Lagrange  multipliers 

A.  (1)  ,  A  (2 )  , A.  (N)  ;  the  new  performance  index  L  can  be  defined 

as  follows: 


L-  "l-x  ^(JV)  Fx(/7)  +  {[x  ’"(A)  *  v  i?j  v(;c)  J  +  A  ^(A-1)  ■  (257  ) 

[<fx(ic)  +Hv{k)  -x(if+l)  ]  +  [Gx(Jc)  *Hv{k)  -x(ic*l)  ]  ’'A  (A+1) } 

It  is  a  well  known  fact  that  minimization  of  the  function  L 
is  equivalent  to  minimization  of  J  when  it  is  subject  to  the 
equality  defined  by  Eqn.  (255) .  In  order  to  minimize  the 
function  L,  one  must  differentiate  L  with  respect  to  each 
component  of  the  vectors  x(k) ,  v(k),  and  A(k)  and  set  the 
results  equal  to  zero.  Thus  we  set 


dL 

dxj  ik) 


-0 


i-l,2,...,b  ;  k’‘l,2,...,N 


dL 

dv^  ik) 


0 


;  k-l,2,...,N-l 


(258) 

(259) 


dL 

dX^ik) 


0 


i>l,2,..  ,  n  ;  k-l,2,...,N 


Dropping  the  subscript  for  simplicity,  Eqns.  (258), 
and  (260)  may  be  obtained  as  follows: 


(260) 
(259) , 


dL 

dxik) 


-0 


($xik)  +(5^A(k+l)  -Xik)  -0  A-1,2,...,  AT-l 


(261) 
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Fx{N)  -  kiN)  -"O  (262) 

R^v(k)  ^  H'^kiki-1)  -0  k~0,l,...,N-l  (263) 

dx(k-l)  ^  Hv(k-l)  -xik)  ^0  k-l,2,...,N  (264) 

After  the  simplification  of  the  equations  just  obtained, 
there  results, 

k(k) -dxik)  ^G'^kik^l)  k~l,2,...,N-l  (265) 

with  the  final  condition 

k{N)^Fx{N)  (266) 

Rearrangement  of  the  terms  of  Eqn  (263)  leads  to  the 
solution  of  v(k) 

vik)  -~R-^H'^X{k*l)  (267) 

The  last  partial  differential  Eqn.  (264)  is  simply  the  state 
scjuation  (see  Eqn.  (255)).  Substitution  of  Eqn.  (267)  into 

Eqn.  (255)  results  in 

x(.k*l)  -Gx{k)  ~  H^kik*!)  ;  x(0)-c  (268) 

In  order  to  obtain  the  solution  to  the  minimization  problem 
we  need  to  solve  Eqns.  (265)  and  (268)  simultaneously. 

Notice  for  the  system  equation,  Eqn.  (268) ,  the  initial 
condition  is  specified,  while  for  the  Lagrange  multiplier 
equation,  Eqn.  (265),  the  final  condition  is  specified.  Thus 
the  problem  here  becomes  a  two-point  boundary  value  problem 


dL 


dx(N) 

dL 

dv(k) 

dL 

■  dk(k} 


-  0 
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(TPBVP) .  If  the  TPBVP  is  solved,  then  the  optimal  values  for 
the  state  vector  and  Lagrange  multiplier  vector  may  be 
determined  and  the  optimal  control  vector,  v(k)  may  be 
obtained  in  the  open-loop  form.  However,  if  one  employs  the 
Ricatti  transformation,  the  optimal  control  vector,  v{k) , 
can  be  obtained  in  the  following  closed-loop  or  feedback 
form: 

v(if)  (269) 

where  is  the  rxn  feedback  matrix. 


Under  the  assumption  that  A{k)  can  be  written  in  the 
following  form 

k{k)  •f(k)x{k)  (270) 

where  P(k)  is  an  nxn  Hermitian  matrix  (see  Eqn.  (250)). 
Substitution  of  Eqn.  (270)  into  Eqn.  (261)  results  in 

^(Jc)x(Jc)  -£?x(k)  +G^j5(Jf+l)  x(k+l)  (271) 

and  substitution  of  Eqn.  (270)  into  Eqn.  (268)  gives 


x{k*l)  -Gxik)  -HR;^H^Pik^l)x(k-^l)  (272) 

Notice  that  Eqns.  (271)  and  (272)  do  not  involve  A.(k)  and 
thus  X(k)  has  been  effectively  eliminated.  The 
transformation  process  employed  here  is  called  the  Ricatti 
transformation.  It  is  of  extreme  importance  in  solving  such 
a  TPBVP.  From  Eqn.  (272) 


IOC 


]  Xik*l)  -Gxik)  (273) 

because  of  the  existence  of  the  inverse  matrix,  Eqn.(273) 
can  be  written  as 


x(Jc+l)  -  lI*HRi^H^P{k*l)  ]  -'^Gxik)  (274) 

By  substituting  Eqn.  (274)  into  Eqn.  (271)  one  obtains 


Pik)  -(?+<5'jP(Jc+1)  +  ] '"G  (275) 

Equation  (275)  is  the  same  as  Eqn.  (250),  it  may  be  modified 
to 

p(k)  -Q*(S'^P(k*l)  (S-G’^P{k->-l)  H[R^*H"P{k-^l)  H]-'^H"P{k*l)  G  (276) 
Equation  (276)  is  called  the  Ricatti  equation.  Referring  to 
Eqn.  (262)  notice  at  k=N 


P(N)x(N)  -K(N)  -Fx(/7) 
P{N)  -F 


(277) 


Hence,  Eqns.  (275)  and  (276)  can  be  solved  uniquely  bac)cward 
from  k=N  to  k=0.  So  one  can  obtain  P(N)  ,  P(N-l),...,  P(0) 
starting  from  P(N)  which  is  known.  By  referring  to  Eqns. 

(270)  and  (274),  the  optimal  control  vector,  v(k)  given  by 
Eqn.  (267)  now  becomes 


v{k)  -  H'^P(k*l)  x(k*l) 

--R^^H'^Pik*!)  [I*HRi^H^P(k^l)]  '^Gxik) 


(278) 
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A  modified  form  of  the  optimal  control  vector,  v(k)  can  be 
given  by 

vik)  H]  --H"P(k*l)  Gxlk)  (279) 

the  same  as  Eqn.  (249) 

5.3.1  -  Evaluation  of  the  Minimum  Performance  Index 

Evaluation  of  the  minimum  performance  index  given  by 
Eqn.  (254)  follows 

dniin"  Itin  j-i-X  '*^'2  52  +  V  ■'■(>:)  v(^)  ]  I  (280) 

Premultiplying  both  sides  of  Eqn.  (271)  by  x^(k)  gives 

x'^'ik)  f(k)  xik)  •x'^(k)  ($x(k)  *x‘^(k)  G  ^i^(A:+l)  x(k*l)  (281) 
Substituting  Eqn.  (273)  and  arranging  terms  leads  to 

X  ’■(k)  f{k)x{k}  -X  ’■(k)  £?x(k)  •►x’'(k*l)  [I*  fl.k*i)  Pl.k*l)  x{k*i)  (282) 

Hence, 

x^(k)  Qxik)  -x^ik)  Px(k)  -x^(k*l)  Pik->-l)  xik*l)  +  (2831 

-X  '^(k*l)  P{k*i)  hr;^^h  ^P{k*i)  x(k*l) 

similarly,  from  Eqn.  (278) 


vik)  ”~r;^ H^Pik*l)  xik^l) 


(284) 


Hence , 


v^(Jc)  R^v(k)-[-x^{k*l)  P{k*l)  h’Ri^]  /?!  "P(Jc+l)x(A  +  l)  j  (285) 

-x"{k-l)  P(k*l)  NR;^H'^f(k*l)  x{k*l) 

By  adding  Eqns.  (283)  and  (285) 
x"^  {k)  Cx{k)  *  v'^ik}  R^v(,k)  -  x  (k)  Px{k)  -  x  ‘^{k*l)  Pik*l)  x  {k*l)  (286) 

By  substituting  Eqn.  (286)  into  Eqn.  (280)  one  obtains 

Jj^  -  -^x^iN)  Fx(N)  ■*■—5^  [x'^ik)  P{k)  x{k)  ~  x ik+1)  Pik+l)  xik*!)  j 
2  2 

-Ix'^iN)  FxiN)  *^[x^{0)P{Q)xiQ)-x"{l)Pil)x{l)^x^il)P{l) 

-X^{2)  Pi2)  x(2)  +  ...  -^-x'iN-l)  P{N-l)x(N-l)  ~  x  (N)  PiN)  X  iN)  ] 
-~x^{N)  Fx(N)  *  ^x'^iO)  P{0)X(0)  -^x^(N)  P{N)xiN) 

(287) 

Notice  from  Eqn,  (277)  P(N)=F.  Hence,  Eqn.  (287)  becomes 

J„in-^X^(0)i5(0)x(0)  (288) 


Thus,  the  minimum  value  of  the  performance  index  J  is  given 
by  Eqn.  (280).  It  is  a  function  of  P(0)  and  the  initial 
state,  x(0) . 


CHAPTER  6 


RESULTS  AND  DISCUSSION 


For  this  thesis  three  different  comparisons  have  been 
analyzed  for  the  nominal  orientation  of  Case  (1).  The  first 
is  a  comparison  of  two  sets  of  actuator  locations  (see  Fig. 
8) ,  this  will  provide  the  most  effective  placement  of  the 
actuators.  The  second  is  a  parametric  study  involving  the 
state  penalty  matrix  and  the  control  penalty  matrix  (see 
Eqns.  (185)  and  (186)),  this  will  yield  the  best  choice  of 
penalty  matrices  under  the  given  conditions.  The  last 
comparison  is  also  a  parametric  study,  through  the  variance 
of  the  sampling  period  (see  Section  5.2)  a  display  of  its 
effect  on  the  system  performance  will  be  produced.  Together 
these  three  comparisons  will  result  in  the  design  of  the 
optimum  cuiitrol  system. 

For  each  comparison  three  different  system 
characteristics  will  be  analyzed.  The  first  characteristic 
is  the  discretized  open  and/or  closed  loop  system  eigenvalue 
and  their  moduli.  In  a  z-domain  the  system’s  eigenvalues  or 
moduli  should  display  the  following  characteristics: 

a)  the  magnitude  should  be  less  than  1  to  ensure  stability 
of  the  system. 

b)  the  magnitude  should  be  as  small  as  possible.  The 
system  with  the  smaller  eigenvalues  usually  displays 
the  better  system  response. 


The  second  analysis  is  the  tabulation  of  the  calculated 
optimum  cost  function  (see  Section  5.3.1).  The  optimum  cost 
function  -also  known  as  the  minimum  performance  index-  is 
calculated  through  the  minimization  of  the  cost  function, 
therefore,  the  obvious  selection  for  an  optim.um  system  is 
the  smallest  possible  value  of  x  x  . 

The  transient  time  response  of  the  rotational  angles, 
modal  amplitudes,  and  the  control  forces  of  each  actuator  is 
the  third  system  characteristic  investigated.  A  transient 
response  refers  to  the  process  generated  in  going  from^  an 
initial  state  to  the  final  state.  Consideration  is  giving 
to  several  elements  of  the  ti^n&ient  response,  they  are: 

a)  the  maximum  o^-ershoot  or  undershoot  which  is  directly 
related  to  the  roouat. cf  the  system.  In  the  case  cf 
the  modal  amplitudes,  the  response  is  due  to  an  initial 
normalized  deflection  of  0.01,  which  corresponds  to  a  1 
meter  deflection.  If  the  response  overshoot  or 
undershoot  is  too  high,  internal  stresses  m.ay  reach  a 
maximum  resulting  in  fracture  or  failure  of  the 
structure.  In  reference  to  the  forces,  a  large  initial 
overshoot  is  associated  with  the  maximum  amplitude  of 
the  control  force 

b)  the  area  under  the  response  curves  is  being 
investigated.  When  comparing  the  rotational  angle 
response  there  is  no  overshoot  or  undershoot.  The 


system  is  subjected  to  an  initial  angular  deflection  of 
0.01  radian.  The  area  under  this  response  and  the 
tangential  envelopes  created  by  the  other  responses  is 
proportional  to  the  control  energy  dissipated; 
therefore,  the  overshoot  and  the  area  should  be  as 
small  as  possible. 

c)  the  rise  time  is  the  time  required  for  the  response 
initially  to  reach  the  equilibrium  value  during  the 
first  cycle  of  the  response. 

d)  the  peak  time  is  the  time  it  takes  the  response  to 
reach  the  first  peak  of  an  overshoot  or  undershoot. 

e)  the  settling  time  is  an  important  characteristic  to 
take  into  consideration.  The  settling  time  is  the  point 
in  time  when  the  response  curve  reaches  and  stays 
within  5%  of  the  ec[uilibrium  value  of  the  system's 
response.  Hence,  one  would  prefer  this  duration  of  time 
to  be  as  short  as  possible. 

f)  the  signal  reconstruction  is  one  other  system 
characteristic  which  is  especially  important  in  the 
comparison  of  the  sampling  periods.  The  transient  time 
response  should  appear  as  smooth  as  possible,  thus 
exhibiting  a  proper  signal  reconstruction. 

In  the  subsequent  sections  of  Chapter  6,  three 
comparative  studies  will  be  discussed. 
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6.1  -  Comparison  of  Actuator  Positions 

For  the  comparison  of  the  actuator  locations,  two  sets 
have  been  chosen.  Set  A  and  A*  (see  Fig,  8) .  The  nominal 
orientation  of  Case  (1),  with  a  sampling  period  of  5 
seconds,  and  the  set  of  penalty  matrices  Sq=0,5  and 
Sr=1.0xl0'’°  is  selected  as  the  standard  system. 

The  system's  closed-loop  eigenvalues  and  optimum  cost 
functions  are  displayed  in  Table  5.  Notice  that  the  minimum 
modulus  and  the  overall  average  modulus  of  Set  A  eigenvalues 
are  smaller  than  those  for  Set  A' .  The  optimum  cost  function 
is  smaller  for  Set  A  than  Set  A'. 

Upon  inspection  of  the  transient  time  responses,  it  is 
seen  that  the  response  for  the  rotational  angles  are  exactly 
the  same  (see  Figs.  11  and  12) .  The  responses  for  the  first 
and  second  modal  amplitude  are  quite  different.  On  the  other 
hand,  the  third  mode's  amplitudal  responses  are  exactly  the 
same  (see  Fig.  13).  Set  A'  first  mode  has  a  higher  overshoot 
and  a  much  larger  settling  time  than  that  of  Set  A. 
Similarly,  the  second  mode  has  a  larger  overshoot, 
undershoot,  curve  area,  rise  time,  peak  time,  and  settling 
time  for  Set  A'  than  for  Set  A.  The  same  follows  for  the 
transient  responses  of  the  forces  for  actuators  1-4  (see 
Figs.  14  and  15) .  The  actuators  5  and  6  are  in  the  same 


position  for  both  sets;  thus,  there  isn't  any  difference  in 
transient  time  response  for  the  two  actuator  forces  (see 
Fig.  16) .  Also,  notice  the  symmetry  of  response  for  the 
actuator  groups  1  &  3,  2  &  4,  and  5  &  6.  This  is  due  to  the 
fact  that  the  actuators  are  placed  in  a  symmetrical  pattern 
around  the  center  of  the  plate  (origin) . 

The  choice  between  actuator  location  sets  is  an  easy 
one.  Obviously,  the  system  of  Set  A  exhibits  the  best 
characteristics-  This  is  an  expected  result  for  a  couple  of 
reasons.  The  first  is:  the  torque  arm  is  a  maximum  for  Set 
A,  thus,  reducing  the  maximum  force  needed  to  control  its 
rigid  rotational  motion.  The  second  is:  Most  of  the  Set  A' 
actuators  are  placed  in  the  vicinity  of  the  nodal  lines  -for 
the  first  three  modes  included  in  the  system  model-  thus, 
possibly  reducing  their  overall  effect. 
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Table  S 

Tbe  Effect  of  Different  Actuator  Locations  on  the 
Closed- Loop  System  Eigenvalues  and  Moduli 


set  A 

.  1 

Set  A> 

Eigenvalues 

Moduli 

Eigenvalues 

Moduli 

Real 

Imag. 

Real 

Imag . 

■B 

0.20871 

0 

0.20871 

0.20872 

0 

0.20872 

2 

0.20873 

. 00177 

0.20873 

0.20873 

.00174 

0.20874 

3 

0.20873 

-.00177 

0.20873 

0.20873 

-.00174 

0.2C874 

■1 

0.22710 

0 

0.22710 

0.22273 

0 

0.22273 

5 

0.25698 

0 

0.25698 

0.26381 

0 

0.26381 

6 

0.79820 

0 

0.79820 

0.79575 

0 

0.79575 

■a 

0.92454 

0 

0.92454 

0.93528 

0 

0.93528 

8 

0.88714 

.36217 

0.9452 

0.89757 

.36118 

0.96752 

9 

0.88714 

-.36217 

0.9452 

0.89757 

-.36118 

0.96752 

10 

0.99443 

0 

0.99443 

0.99443 

0 

0.99443 

11 

0.99443 

0 

0.99443 

0.99443 

0 

0.99443 

12 

0.99443 

0 

0.99443 

0.99443 

0 

0.99443 

XjP„X  =0.0086919 

.a  0  D 

bbbebebbsbi 

-  Optimum  Cost  Function 
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Figure  11.  Transient  Time  Response  Of  the  Rotational 
Angles*  for  Comparison  of  Actuator  Positioning,  Set  A 
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Figure  13.  Transient  Time  Response  Of  the  Modal  Amplitudes' 
for  Comparison  of  Actuator  Positioning,  Set  A  &  A' 
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Figure  14,  Transient  Time  Response  of  the  Actuator  Force  At 
162  for  Comparison  of  Actuator  Positioning,  Set  A  6  A* 
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Figure  15.  Transient  Time  Response  of  the  Actuator  Force  At 
3  i  4  for  Comparison  of  Actuator  Positioning,  Set  ASA' 
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Figure  16.  Transient  Time  Response  of  the  Actuator  Force  At 
566  for  Comparison  of  Actuator  Positioning,  Set  ASA' 
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6.2  -  Comparison  of  Penalty  Matrices 

For  the  parametric  study  of  the  penalty  matrices, 
several  system  variables  were  kept  constant;  they  were 
chosen  as:  the  nominal  orientation  of  Case  (1),  Set  A  for 
the  actuator  locations,  and  a  sampling  period  of  5.0 
seconds.  The  variance  of  the  system  penalty  matrices  is 
handled  by  changing  the  control  penalty  matrix  while  holding 
the  state  penalty  matrix  constant,  and  vice-versa.  The 
penalty  matrices  are  defined  by  Eqns.  (185)  and  (186)  where, 
and  a^^Sr.  The  sets  of  penalty  matrices  are: 

a)  Set  1  -  Sq,=0.5,  Sr, =1 . 0x10’’° 

b)  Set  2  -  Sq,=0.5,  Srj^l .  OxlO’’^ 

c)  Set  3  -  Sqj^O.OOl,  Sr, =1.0x10*''° 

d)  Set  4  -  Sqj^O.OOl,  Sr,=l .  OxlO*''^ 

Note  that  in  this  study  only  block  diagonal  penalty  matrices 
-where  all  the  elements  are  the  same-  are  considered. 

The  eigenvalues,  moduli,  and  optimum  cost  function  for 
each  set  are  displayed  in  Table  6.  The  eigenvalues  and 
moduli  of  Sets  1,  2,  and  4  are  very  close,  whereas,  the 
first  modulus  of  Set  3's  eigenvalues  is  noticeably  larger. 
This  helps  to  narrow  the  system  comparison  down  to  the  three 
sets  1,  2,  and  4.  After  taking  a  second  look  at  their 
eigenvalues,  it  is  noticed  that  the  moduli  of  Set  4's 
eigenvalues  are  slightly  larger  than  those  of  the  two  other 


sets;  but  the  optimum  cost  function  for  Set  4  is  smaller 
than  those  for  Sets  1  and  2. 

For  a  complete  comparison  the  transient  response  is 
viewed.  The  responses  for  the  system’s  rotational  angles 
pictured  in  Figure  17  show  no  appreciable  difference 
regardless  of  the  penalty  matrices  used.  In  addition,  there 
is  no  significant  difference  between  the  yaw,  pitch,  and 
roll  responses.  The  transient  response  of  the  modal 
amplitudes  show  a  slight  variation  as  a  function  of  the 
penalty  matrices  used.  For  the  first  and  second  mode  there 
is  virtually  no  difference  in  responses  for  all  combinations 
of  penalty  matrices  (see  Figs.  18  and  19) .  The  third  mode 
responses  display  a  deviation  towards  a  better  response  for 
Set  4  (see  Fig.  20) .  Similarly,  all  the  actuator  force 
responses  (see  Figs.  21~26)  exhibit  better  characteristics 
for  the  case  of  Set  4  penalty  matrices. 

The  determination  of  the  best  combination  of  the  state 
penalty  matrix  and  control  matrix  is  a  little  more  difficult 
than  the  previous  choice  of  positioning.  After  taking  into 
consideration  the  system  characteristics,  Set  4  (Sq=0.001, 
and  Sr=l .  0x10"’^)  was  chosen.  Even  though  the  Set  4  system 
eigenvalues  are  not  the  smallest  ones  considered,  its 
optimum  cost  function  is  the  least  and  its  response  is  the 
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Table  6 

Effect  of  Different  Penalty  Matrix  Cciabinations  on  the 
Closed-Loop  System  Eigenvalues  and  Moduli 
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Figure  17 .  Transient  Time  Response  of  the 
Angles'  for  Comparison  of  Penalty  Matrices 
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Figure  21.  Transient  Time  Response  of  the  Actuator  Force  At 
1  for  Comparison  of  Penalty  Matrices 
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Figure  22.  Transient  Time  Response  of  the  Actuator  Force  At 
2  for  Comparison  of  Penalty  Matrices 
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Figure  24.  Transient  Time  Response  of  the  Actuator  Force  At 
4  for  Comparison  of  Penalty  Matrices 
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For  the  comparison  of  sampling  periods,  several  system 
characteristics  kept  constant  were  chosen  to  be:  the  nominal 
orientation  of  Case  (1) ,  the  actuator  placement  of  Set  A, 
and  the  penalty  matrices  of  Set  1.  The  different  sampling 
periods  chosen  for  the  comparison  were  2.5,  5.0,  10.0,  and 
40.0  seconds,  respectively.  The  latter  was  an  approximate 
value  chosen  from  Table  4  (see  39.95  seconds)  as  a 
representation  of  system  degradation. 

Each  set's  discretized  open  and  closed  loop  system 
eigenvalues  and  moduli  and  their  optimum  cost  functions  are 
shown  in  Table  7.  The  open-loop  eigenvalues  of  all  four 
sampling  periods  show  signs  of  instability  by  containing 
moduli  greater  than  1.  This  is  a  confirmation  that  the 
uncontrolled  (local  horizontal)  orientation  of  Case  (1)  is 
unstable.  When  the  closed-loop  eigenvalues  of  each  sampling 
period  are  examined,  an  unexpected  result  surfaces.  Almost 
all  the  moduli  of  the  40.0  second  eigenvalues  are  slightly 
smaller  than  those  for  the  other  sampling  periods.  This 
results  from  the  fact  that  the  calculated  unacceptable 
sampling  period  is  39.95  seconds  and  not  40.0  seconds,  i.e., 
roundoff  error.  The  sampling  period  of  40.0  seconds  was 
employed  to  provide  an  equivalent  final  time  period  with  the 
ocher  sampling  periods,  in  other  words,  2.5,  5.0,  10.0  are 
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multiplies  of  40.0.  The  moduli  of  the  last  eigenvalue  of 
Set  IV  is  very  close  to  the  limit  of  1,  thus  exhibiting  a 
trend  towards  degradation.  Comparison  of  the  optimum  cost 
function  indicates  a  minimum  for  the  2.5  second  sampling 
period. 

The  transient  time  response  for  the  rotation  angles  are 
the  same  for  all  the  sampling  periods  (see  Fig.  27) .  When 
the  transient  response  of  the  modal  amplitude  are  considered 
the  first  mode's  maximum  overshoot  decreases,  its  undershoot 
increases,  as  the  sampling  period  increases  (see  Figs.  28- 
31) .  The  same  relationship  holds  for  the  response  of  the 
second  mode.  The  third  mode's  response  displays  an 
undershoot  for  the  10.0  and  40.0  seconds  sampling  periods,  a 
characteristic  not  shown  before  in  this  study.  The  entire 
modal  responses  for  the  40.0  second  sampling  period 
represent  a  degradation  in  comparison  with  the  modal 
responses  for  the  other  sampling  periods.  For  the  case  of 
the  transient  force  responses,  the  maximum  overshoot 
decreases  as  the  sampling  period  increases  (see  Figs.  BI¬ 
SS).  For  the  40.0  second  case  there  is  no  overshoot  but  a 
large  undershoot  and  a  very  choppy  response,  again 
associated  with  system  performance  degradation. 

Initially,  the  selection  of  the  sampling  period  can  be 
reduced  to  just  three  periods,  2.5,  5.0,  and  10.0  seconds. 
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for  the  obvious  reason  that  the  sampling  period  of  40 
seconds  is  associated  with  performance  degradation.  However, 
it  is  not  an  easy  task  to  decide  on  the  best  sampling  period 
from  the  three  other  sampling  periods.  Effectively,  2.5 
seconds  can  be  released  from  consideration  due  to  the  larger 
moduli  of  the  eigenvalues  and  overshoot.  The  decision 
between  the  last  two  sampling  periods  is  a  more  difficult 
one.  The  sampling  period  of  5.0  seconds  has  a  higher 
overshoot,  a  smaller  undershoot,  a  smaller  optimum  cost 
function,  and  a  smooth  response  curve.  The  sampling  period 
of  10.0  seconds  displays  a  lower  overshoot,  a  larger 
undershoot,  smaller  eigenvalues,  and  a  choppy  response. 
Finally,  taken  into  consideration  that  signal  reconstruction 
is  very  important,  the  sampling  period  of  5.0  seconds  can  be 


chosen. 


Table  7 

The  Effect  of  Different  Saapling  Periods  on  the  Open 


Table  7 
continued 
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Figure  27.  Transient  Time  Response  of  the  Rotational 
Angles’  for  Comparison  of  Sampling  Period,  5.0  seconds 
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Figure  29.  Transient  Time  Response  of  the  Modal  Amplitudes 
for  Comparison  of  Saunpling  Period,  5.0  seconds 
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Figure  32.  Transient  Time  Response  of  the  Actuator  Forces 
At  All  Positions  for  Comparison  of  Sampling  Period,  2.5 
seconds 
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CHAPTER  7  -  CONCLUSIONS 

A  mathemdi.  leal  model  to  predict  the  dynamics  of  a 
flexible  orbiting  platform  was  developed.  Under  the 
assvimption  that  the  linear  system  is  completely  observable, 
the  optimal  control  laws  are  developed  for  the  case  where 
the  observational  data  is  collected  on  a  sampled  basis, 
i.e.,  a  discrete  time  data  system. 

Attitude  and  shape  control  of  the  platform  was  assumed 
to  be  provided  by  the  placement  of  point  thrust  actuators 
perpendicular  to  the  main  surface  and  the  edge  of  the  plate. 
Their  effects  on  the  system's  motion  were  modeled  to  the 
first  order.  Controllability  for  the  system  was  verified  for 
two  sets  of  actuator  locations.  An  application  of  the  linear 
quadratic  regulator  (LQR)  techniqfue  in  a  discrete-time 
domain  yields  the  optimum  control  law  feedback  gains. 

A  comparison  of  the  performance  of  the  different  sets 
of  actuator  locations  resulted  in  the  best  choice  of 
actuator  positioning.  Two  parametric  studies  were  conducted 
to  show  the  effect  of  varying  the  state  penalty  matrix  and 
the  control  penalty  matrix,  and  the  effect  of  changing  the 
sampling  period  on  the  transient  performance  of  the  system. 

Generally,  when  comparing  all  of  the  system 
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characteristics,  one  must  search  for  the  system  with  the 
maximum  number  of  positive  characteristics.  Unfortunately, 
it  is  not  always  a  clear  cut  decision  as  to  which 
combination  of  characteristics  designate  the  best  system. 
Often  a  system  may  exhibit  some  value  tradeoffs,  it  may 
possess  several  minimum  values  along  with  several  maximum 
values.  In  this  thesis  particular  emphasis  is  placed  on  the 
quality  of  the  transient  response.  With  this  in  mind  the 
following  conclusions  have  been  made. 

When  deciding  upon  the  best  placement  of  the  actuators 
on  the  main  surface,  two  items  should  be  kept  in  mind.  The 
first  is  that  the  actuators  should  be  placed  such  that  there 
is  the  maximum  distance  between  the  origin  and  the  actuator 
locations.  This  creates  a  maximum  torque  arm  which  reduces 
the  maximum  force  needed  to  control  the  system’s  rigid  body 
motions.  The  second  is  that  the  actuators  should  be  placed 
as  far  away  as  possible  from  the  nodal  lines  of  the 
fundamental  and  lower  frequency  elastic  modes  (with  which 
most  of  the  elastic  energy  is  associated) . 

For  the  determination  of  the  best  choice  of  the  penalty 
matrices,  one  usually  chooses  a  large  state  penalty  matrix 
and  a  small  control  penalty  matrix  within  the  limits  of 
control  saturation  levels.  These  matrices  should  be  chosen 
such  that  they  optimize  the  performance  index  and  produce 
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the  best  transient  time  response. 

Since  the  sampling  period  is  directly  related  to 
discretization  and  the  signal  reconstruction  of  the  system, 
careful  consideration  should  be  given  to  its  choice. 
Initially,  when  deciding  upon  a  sampling  period,  one  should 
take  into  consideration  the  sampling  period  theorem  which 
involves  the  use  of  the  open-loop  system  eigenvalues  to 
calculate  the  unacceptable  sampling  periods  for  which  the 
system  is  not  controllable.  Any  sampling  period  values  near 
the  values  of  the  unacceptable  periods  should  not  be 
utilized,  for  they  can  lead  to  system  performance 
degradation.  The  sampling  period  should  be  chosen  such  that 
the  sampling  period  is  not  too  small  to  avoid  excessive 
handling  and  accumulation  of  data  for  the  onboard  computer 
system.  On  the  other  hand,  the  sampling  period  should  not  be 
too  large,  otherwise  the  system’s  transient  time  response 
will  appear  choppy,  thus,  failing  to  provide  a  good  example 
of  signal  reconstruction. 


SUGGESTIONS 


CHAPTER  8  - 

In  terms  of  new  actuator  positioning,  different 
locations  for  the  actuators  placed  on  the  edge  of  the  plate 
should  be  investigated.  In  addition,  the  number  of  actuators 
could  be  varied,  in  an  effort  to  improve  system  performance. 
There  are  some  system  tradeoffs  to  consider  in  the  event 
that  the  number  of  actuators  is  increased  or  decreased.  If 
the  number  of  actuators  are  decreased  there  may  be  some  loss 
in  system  controllability  and  robustness.  On  the  other  hand, 
if  the  number  of  actuators  are  increased,  the  system's  total 
mass  will  then  be  increased,  thus,  creating  an  increase  in 
the  system's  operational  costs. 

One  other  possible  improvement  to  actuator  modeling 
would  be  incorporation  of  lumped  masses  to  represent  the 
actuator  mass  at  its  particular  location.  Implementation  of 
GTSTRUDL  or  similar  finite  element  methods  for  the 
recalculation  of  modal  frequency  and  deflections  could  than 
be  performed. 

When  considering  possible  improvements  in  the  choice  of 
penalty  matrices,  a  more  indepth  comparison  should  be 
performed.  The  split  weighting  penalty  matrix  approach  could 
be  applied,  with  a  possible  variance  of  each  individual 
matrix  element.  This  could  result  in  a  smaller  optimum  cost 


function. 


For  this  thesis  the  system  considered  was  assumed  to  be 
completely  observable  and  deterministic,  it  is  an  inadequate 
assumption.  To  improve  system  modeling  the  observer  matrix 
should  be  developed  and  the  system  should  be  considered 
stochastic  (i.e.,  experience  external  disturbances  and/or 
noise) .  The  linear  quadratic  Gaussian  (LQG)  technique  can  be 
applied  to  provide  the  control  law  and  estimator  for  a 
stochastic  system  with  a  Kalman  filter  to  act  as  a  screen. 

One  last  improvement  to  the  system  corresponds  to  the 
initial  assumptions  about  the  material  properties  of 
composite  graphite.  For  this  thesis  the  material  is  assumed 
to  be  isotropic.  This  is  an  questionable  assumption  for 
reinforced  composite  graphite.  Generally,  the  material  is 
composed  of  graphite  fibers  bonded  by  a  resin  epoxy.  Even  if 
the  fibers  were  aligned  with  a  9o;  0°  or  4  5'  orientation, 

the  material  at  best  can  be  assumed  to  be  orthotropic . 

The  system’s  various  material  constants  would  then  need  to 
be  redefined.  These  material  properties  affect  the  stress- 
strain  relationship  of  Eqn.  (125)  and  create  cross  terms  in 
the  plate  vibrational  equation,  thus  affecting  the  vaiues  of 
the  natural  frequency.  The  new  values  for  Young's  Modulus, 
Poisson's  ratio,  and/or  Shear  Modulus,  (i.e.,  E^,  E^,  E^, 
etc.)  could  then  be  submitted  to  GTSTRUDL  for  recalculation 
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of  the  modal  frecpaencies  and  deflections. 

The  results  of  the  present  study  could  then  be  compared 
to  the  previous  LQG  results  with  the  thin  flat  plate  and  the 
shallow  spherical  shell  for  corresponding  actuator 
positions. 
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Jc^^-Sa^  +  h^Ci-v)  -5.3x10* 
i:,g-15ar-3jb(l-v)  -1,290 
ic77-20a^  +  4i?^  (1-v)  -2.28x10® 
irgj-lOh^-a^d-v)  -9 .3x10* 
icg3-3  0t>/r-3a(l-v)  --3,210 
iCgg-Sjb^  +  a^d-v)  -5.7x10* 
Jc0s-i5i>/r-3ad-v)  -1,290 
lC8,-15ab-l  .5x10® 
i:e8"20b^  +  4a*  (l-v)  -2.28x10® 
iCji  - -l5ar  +  l5bv  -  3bd-v)  --840 
Jr52  --3  0b/r-3ad-v)  --3,210 
Jc93-60/rd30r2-30v-42  (l-v)  -51.6 
JCg^- -15a/ir  +  3bd-v)  --1,290 
icgg  - -I5jb/r- 15av  - -1 , 29  0 
iCjg  -  -30/r^  -  30r^  +  30v  *42  d-v)  -  -21 . 6 
iCg, --30ar-l5bv-3bd-v)  --3,660 
Jc9g--30i)/r-15av-3ad-v)  --3,660 
iCgg  -  15b/r- I5av  - 3a  (l-v)  -840 


+  (1-v)  -S.lxlO* 

JCio_3--15ar  +  3b(l-v)  -1,290 
iCj^o  4- 10a^-4jb2  (i-v)  -7  .2x10* 
k^g_f-l5ar-l5bv  -3b{l-v)  -840 
ic^o  .,-10a2-i)2(l-v)  -9.3xi0‘ 
^io,9"-30ar-3i?(l-v)  --3,210 
^io,io*20a2  +  4i32  (i-v)  -2.28x10^ 

•^10,2  ”  -^10,5  "  -^10,8  "  -^11,1  "  4  “  -^11,7  “ 

iCii  j-Si^  +  a^d-v)  -S.lxlO* 
^33_3-15jb/r2-3a(l-v)  -1,290 
Jc^3_5  -  lOjb^  -  a‘  (1-v)  -  9  -iXlO* 

-  30jb/r  +  I5av  +  3a  (1-v )  -3,660 
iCj^3_a-10i?2-4adl-v)  ”1 .2x10* 

g - -I5b/r  +  15av  +  3a (1-v)  --840 
^11.10“  “15ai?v- -1.5x10* 
^u,ii“20i?d4adl-v)  -2.28x10* 
^i2.i“15a/r-3l?(l-v)  -1,290 
^12.2“  -15i?/r  4-3a  (1-v)  --1,260 
^i2.3"~30/r2-30r2  +  3  0v-K42  (1-v)  -  -21.6 
^12.4  "  ^5ar- I5i?v  +  3a  (1-v )  -840 
ic32  5- -30ib/r-3a(l-v)  --3,210 
i:i2,«-6  0/rd30r2-30v-42  (1-v)  -51.6 
ic^2  -,-30ar  +  3ii(l-v)  -3,210 
Jc.j  g- -15l?/r  + 15av -3a(l-v)  --840 
'^i2,9“30/r2-30r2-30v-42  (1-v)  --38.4 
iCjj.io  -  30ar  +  I5fc>v  +  3i?(l -V )  -3,660 
-^12,11 "  -30i7/r-15av  -  3a(l-v)  --3,660 
JCi,  i,-15i)/r- 15av-3a(l-v)  -840 
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Appendix  B  ~  GTSTRUDL  Program  and  a  Summary  of  the 

Implemented  Commands  for  the  Calculation  of 
Modal  Frequencies  and  Shape  Patterns 

Because  the  model  being  considered  is  a  free-free 
structure,  GTSTRUDL  recognizes  the  generated  stiffness 
matrix  as  a  representation  of  an  unstable  structure. 
Consequently,  it  will  not  formulate  a  complete  dynamic  modal 
analysis.  For  a  calculation  of  the  modal  frequency  values, 
the  command,  'COMPUTE  RIGID  BODY  MODES',  must  be 
implemented.  Since  the  rigid  body  modes  must  be  computed, 
this  eliminates  using  the  eigenproblem  solver,  'SUBSPACE 
ITERATION',  for  it  can  not  calculate  those  modes.  There  are 
two  other  eigenproblem  solution  methods  available: 
'TRIADIAGONALIZATION'  and  'GTLANCZOS'.  Triadiagonal ization 
can  not  make  the  correct  calculations,  because  it 
incorporates  the  Inverse  Iteration  method,  which  inverts  the 
matrix  resulting  in  problems  with  singularity.  Therefore, 
for  the  solution  of  the  eigenvalues  the  GTLANCZOS  method 
must  be  used;  it  uses  the  Lanczos  eigenvalue  solver  method. 

The  eigenvalues  are  obtained  by  using  two  GTSTRUDL 
commands,  'PERFORM  EIGEN  ANALYSIS'  and  'LIST  EIGENVALUES'. 
The  final  values  for  the  frequencies  are  approximated 
graphically.  By  using  the  'LIST  EIGENVECTORS'  command  the 
eigenvectors  can  be  computed.  This  gives  the  normalized 
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deflections  at  each  node.  The  modal  shape  patterns  are 
acquired  through  the  combination  of  the  'LIST  EIGENVECTORS' 
and  'PLOT  PLANE  MODE  SHAPE  #XZ  PROJECTION'  commands.  The 
second  command  plots  a  diagram  of  the  normalized  deflections 
of  the  XZ  plane.  The  YZ  plane  was  another  plane  orientation 
viewed;  for  this  orientation  the  Z  axis  is  perpendicular  to 
the  major  surface.  Table  3  represents  the  modal  shape 
patterns  derived  from  GTSTRUDL. 

In  addition  to  the  previous  commands,  there  are  some 
others  that  check  the  performance  criteria  and  calculation 
values:  'ORTHOGONALITY',  ' STURM  SEQUENCE ' ,  and  'ERROR 
ESTIMATE' .  The  ORTHOGONALITY  check  has  been  previously 
mentioned  in  Section  3.4.5.  The  STURM  SEQUENCE  check 
determines  whether  the  correct  number  of  modes  are  within  a 
specified  range.  This  calculation  requires  the  decomposition 
of  a  matrix  that  has  an  order  and  banding  equal  to  the 
system  stiffness  matrix.  Si.nce  the  stiffness  matrix  for  this 
system  is  very  large,  it  was  assumed  this  check  would  need 
too  much  memory  space  and  could  not  be  afforded.  Similarly, 
was  the  situation  with  the  ERROR  ESTIMATE  process.  This 
command  computes  an  error  for  each  mode  as  follows: 

\[K]  {4^,} 

![/:]  {<5,}! 


The  ERROR  ESTIMATE  command  compares  this  calculate  value, 
Cj,  to  the  set  'TOLERANCE'  value  -default  of  1.0x10*-  to 
determine  whether  the  percentage  of  error  is  within  a 
reasonable  range.  For  most  situations  the  default  value  i 
appropriate,  but  may  be  changed;  this  may  have  a  major 
effect  on  solution  time. 
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Appendix  C  -  The  System  Matrices  for  the  Nominal 

Orientations  of  Case  (2^  and  Case  (3^ 


Case  (2)  -  Ix=Iv=Iz/2 


State  Matrix,  [A] : 


U],- 


0  0 


ex6 

0 

0 

-3 

0 

0 

0 


[I] 


0 

0 

0 

-2402 

0 

0 


0 

0 

0 

0 

-49  i9 
0 


0 

0 

0 

0 

0 

-7667 


0 

-2 

0 


6x6 


1 

0 


(C 


Control  Effort  Matrix,  [B]{U): 


[sj  {U] 


2,  A 


[BJ  {U} 


2.  A' 


.  1605 

-.321 

-.1605 

.  321 

0 

.  321 

.  1605 

-.321  - 

.  1605 

0 

0 

0 

0 

0 

.  1605 

- . 2221 

.  2221 

-.2221  . 

2221 

0 

-  .  189 

.  189 

-  .  189 

.  189 

0 

.  1842 

.  1842 

.1842 

1842 

0 

.  1605 

-.214 

- . 1605 

.214 

0 

.214 

.1605 

-.214  - 

- .  1605 

0 

0 

0 

0 

0 

.  1605 

- . 1545 

.  1545 

- . 1545 

.  1545 

0 

-  .  05556 

.  05556 

-.05556 

.  05556 

0 

.  0306 

.  0306 

.  0306 

.0306 

0 

0 

0 

1605 

0 

0 

0 


'A, 

A. 


0 

0 

L6[ 

0 

0 

0 


A, 


.1) 


MC.2) 


(C.3) 
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State  Matrix,  [A]: 


[A],-  0 


0 

0 

0 

0 

0 

0 


0 

0 

0 

-2401 

0 

0 


0 

0 

0 

0 

-4948 . 5 
0 


0  0 
0 

0  0 
0 

0  ; 
-7666  0 


.  0 

:  (C.4) 


Control  Effort  Matrix,  [B]{U}; 


IB]  {m3,A- . 


.  1605 

-.321 

- . 1605 

.  321 

0 

0 

0 

0 

0 

.  1605 

-.321 

-.1605 

.321 

.1605 

0 

-.2221 

.2221 

-.2221 

.2221 

0 

-.189 

.189 

-.189 

.  189 

0 

.1842 

.1842 

.  1842 

.  1842 

0 

A.  (C.5) 

< 


[B] 


.  1605 

-.214 

-.1605 

.214 

0 

0 

0 

0 

0 

.  1605 

-.214 

-.1605 

.214 

.  1605 

0 

-.1545 

.  1545 

-.1545 

.  1545 

0 

-.05556 

.05556 

-.05556 

.05556 

0 

.  0306 

.  0306 

.0306 

.0306 

0 

(C.6) 
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